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SUMMARY 


The electromagnetic modeling of packages and interconnects plays a very important 
role in the design of high-speed digital circuits, and is most efficiently performed by 
using computer-aided design algorithms. In recent years, packaging has become a 
critical area in the design of high-speed communication systems and fast computers, and 
the importance of the software support for their development has increased accordingly. 
Throughout this project, our efforts have focused on the development of modeling and 
simulation techniques and algorithms that permit the fast computation of the electrical 
parameters of interconnects and the efficient simulation of their electrical performance. 



PROJECT DESCRIPTION 


The development of efficient and accurate computer-aided design tools is essential 
for the implementation of high-speed digital circuits used in computer systems and 
communication networks. With current trends in which network complexity and signal 
speed keep increasing, problems associated with signal integrity such as crosstalk, 
distortion, losses can compromise the overall electrical performance of computers and 
communication systems. 

Presently, industrial needs for computer support in network design is increasing 
rapidly; however the availability of design and analysis tools capable of handling the 
complexity and volume of manufactured systems lags seriously. Future winners in the 
competitive world of high-speed communications will have to possess sophisticated 
packaging analysis and design tools and system performance will be more and more 
determined by the adequacy of these tools. 

Because of their important role in the design process, CAD tools must offer certain 
essential features such as speed, accuracy, availability of extensive library and good user 
interface. CAD tools will insure minimum manufacturing cost, faster turnaround time 
and more reliable hardware for production. To achieve the required performance, 
optimization routines must be made available and efficient algorithms must be 
implemented to guarantee speed and good user interface, friendliness and portability of 
the software. 

Our group at the University of Illinois has been involved during the past three years 
in the development of modeling and simulation techniques for interconnects and 
packages. In the parent project, the main emphasis of our work is directed toward the 
modeling of complex interconnect structures as well as the simulation of interconnects. 

The development of efficient and accurate computer-aided design tools is essential 
for the implementation of high-speed digital circuits used in computer systems and 
communication networks. With increasing clock rates and reduced circuit sizes, 
electromagnetic phenomena such as crosstalk noise, distortion, ground bounce will 
become more pronounced in both circuit boards and chip environment. These effects 
seriously degrade the signal integrity of high-speed networks and compromise the 
overall network performance. 



Development of Interconnect Modeling Techniques 

In many situations, one factor that contributes to the increased computation time in 
the calculation of the electrical parameters of the Green's function in the spatial domain, 
which represents the vector field produced by an infinitesimal dipole placed over the 
dielectric substrate layer backed by a ground plane. An efficient method to compute the 
2-D and 3-D capacitance matrix of multiconductor interconnects in a multilayered 
dielectric medium was developed in our group. The method is applicable to conductors 
of arbitrary polygon shape embedded in a multilayered dielectric medium with possible 
ground planes on the top or bottom of the dielectric layers [2], [7]. The method has been 
extended to the computation of equivalent capacitance of via structures in multilayer 
environment [5]. 

In the time domain analysis, the ability to model fine features, e.g., wire bonds, is an 
important requirement and is unavailable in the conventional finite-difference time- 
domain (FDTD) approach unless a very high density of discretization is employed. The 
FDTD method is one of the most widely used techniques. In the FDTD, the derivative 
operators are replaced with the central difference operators, which preserves the second 
order accuracy. Hence, its extension to the general nonuniform grids is not possible 
without losing the second order accuracy or reformulating in terms of the curvilinear 
coordinates. However, in solving any practical problems nonuniform grids are highly 
desirable due to the limitation of computer resources. 

An efficient way to implement the surface impedance boundary condition (SIBC) for 
the finite-difference time-domain FDTD method was introduced [4], Surface impedance 
boundary conditions are first formulated for a lossy dielectric half-space in the frequency 
domain. The impedance function of a lossy medium is approximated with a series of 
first-order rational functions. Then the resulting time-domain convolution integrals are 
computed using recursive formulas which are obtained by assuming that the fields are 
piecewise linear in time. Thus the recursive formulas derived are second-order accurate. 
The preprocessing time is eliminated by performing a rational approximation on the 
normalized frequency-domain impedance. This approximation is independent of material 
properties. 


Simulation of Interconnects 

In the real world of electronic packaging, transmission lines are more likely to be 
nonuniform and may include discontinuities such as bends, tapers and transitions, hence, 
the standard simulation tools for uniform lines can no longer be used to analyze them. 
Presently, a number of methods are available for the simulation of coupled transmission 
lines that are used to model interconnects. In the past two years, we have carried out a 
systematic comparison of these methods with a view to developing an approach that 
would be optimal in terms of both accuracy and efficiency. This has led to the 
development of a transient simulation method based on the difference approximation 
which has the highly desirable feature that it can be conveniently incorporated in a 
circuit simulator [6]-[7], [10]-[13]. This approach not only outperforms the standard 
scattering parameter method, but is very accurate and computationally efficient as well. 
Software designers at Cadence Design Systems and Intel have recently implemented this 
method in the latest circuit simulators. 

The problem of distributed line simulation was analyzed, and the optimal method, 
which results in the maximum efficiency, accuracy and practical applicability, was 
developed. The method is applicable to transmission lines characterized by frequency or 
time-domain data samples. The resulting line model can be directly used in a circuit 
simulator. The efficiency of the optimal method allows for the accurate transient 
simulation of real circuits containing thousands of lossy coupled frequency-dependent 
nonuniform lines surrounded by nonlinear active devices with virtually no increase in the 
simulation time compared to that for the simple replacement of interconnects with 
lumped resistors. 

As components of the optimal method, the following novel techniques were 
introduced: 

- The system model for uniform and nonuniform lines which simplifies analysis of 
distributed networks; the open-loop device model for uniform and nonuniform lines 
which relates voltages and currents at the line terminals via the simplest possible transfer 
functions and time-domain responses; 

- The indirect numerical integration-a class of numerical integration methods, which has 
ideal accuracy, convergence and stability properties, 



- The difference approximation, a general method for applying numerical integration to 
systems characterized by discrete data samples was developed and put in a matrix form. 

- The matrix delay separation from the matrix propagation function, which avoids the use 
of frequency-dependent modal transformation matrices; 

- The relaxation interpolation method, which allows for an accurate and efficient 
approximation of line responses in the time and frequency domains, automatically 
reduces the approximation order depending on the original function and eliminates 
spurious positive poles. 

The complete set of frequency-domain relationships between the matrix Z, Y and S 
parameters were derived and the direct interpolation-based complex rational 
approximation method for transient simulation of macromodels for complicated multiport 
interconnects (such as IC packages and connectors) was developed. The direct 
interpolation-based method was applied to the automatic generation of lumped 
equivalent-circuit models of multiport EM systems (with Dan). 

Approximation Techniques for Circuit Analysis 

Our research has also been focused in developing a unified methodology of model- 
order reduction techniques for circuit and interconnects simulation. The following three 
classes of model-order reduction methods: moment-matching technique, Krylov subspace 
techniques, and reduced optimum approximation have been studied and their applications 
for efficient circuit simulations have been identified. 

The moment-matching technique has been shown to be very effective for generating 
low order models for linear lumped and distributed systems. The method is useful for 
systems whose main features can be retained by the first few orders of reduced system 
models. These include the response estimations of linear lumped networks of medium 
complexity, wave propagation in transmission lines with short delays and diffusion 
process in p-n junctions [20]. 

Krylov subspace based methods such as Amoldi algorithm and Lanczos algorithm are 
effective and robust in generating reduced-order model of large complex systems 
described by ordinary differential or difference equations. The methods are important in 
obtaining reduced-order models to systems that can be characterized by relatively higher 



order models. The methods avoid the construction of ill-conditioned moment matrices 
and the loss of information contained by the eigenvalues of the systems with smaller 
magnitudes when higher models are sought. The methods are suitable for analyzing large, 
complex lumped networks. 

An optimum approximation in conjuction with model order reduction techniques such 
as balanced representation and aggregation methods represents a very effective 
methodology for generating reduced-order models of complex interconnects (infinite 
dimension systems). The approximation method is applied to the distributed systems to 
drive high-order finite models as an intermediate stage and then using balanced 
transformation and aggregation method lower-order models are generated. 
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Optimal Transient Simulation of Transmission Lines 

Dmitri Borisovich Kuznetsov, Student Member, IEEE, and Jose E. Schutt-Aine, Member, IEEE 


Abstract — This paper presents an attempt to formulate a 
high-level description of the optimal transmission line simulation 
method. To formulate the optimal approach, most significant 
aspects of the problem are identified, and alternative approaches 
in each of the aspects are analyzed and compared to find the 
combination that results in the maximum efficiency, accuracy 
and applicability for the transient analysis of digital circuits. 
The practical implementation of the optimal method for uniform 
multiconductor lossy frequency-dependent lines characterized by 
samples of their responses is outlined. It Is shown on an extensive 
set of runtime data that, based on the optimal approach, the 
accurate line modeling in a circuit simulator is as efficient as the 
simple replacement of interconnects with lumped resistors. 

I. Introduction 

T he PROBLEM of transmission line simulation gained 
special importance with the development of high-speed 
digital electronics. As transient times become faster, the trans- 
mission line behavior of electronic interconnects starts to 
significantly affect transient waveforms, and accurate mod- 
eling of on-board and even on-chip interconnects becomes 
an essential part of the design process. The complexity of 
contemporary digital circuits necessitates the simultaneous 
simulation of thousands of lossy coupled frequency-dependent 
lines surrounded by thousands of nonlinear active devices. 
Lines to be simulated may be characterized by measured or 
electromagnetically simulated samples of their responses. 

There are two approaches to the interconnect simulation. 
The first approach creates macromodels for linear subcircuits 
that may contain many transmission lines and other linear 
elements [1], [2]. This paper discusses the second approach, 
in which each multiconductor line is treated as an individual 
circuit element. 

The problem of the line simulation involves several areas of 
science, such as electromagnetics, computational mathematics, 
and circuit and system theories. The solution of the problem 
is straightforward in the sense that all of the components 
involved are well known and only have to be combined 
together. The integration of areas, however, is a difficulty 
that keeps the problem open and accounts for the diversity 
of developed methods. 

This paper presents an original attempt to identify the 
components of the problem and to formulate a high-level de- 
scription of the method that provides the maximum efficiency, 
accuracy and applicability for the transient analysis of digital 
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circuits. Such an approach allows one to more accurately 
assess and compare the performance of numerous and diverse 
line simulation methods. 

The next section presents the formulation of the optimal 
approach and discusses many of the existing line simulation 
methods. Section III outlines the authors implementation of 
the optimal method and presents numerical verification of the 
method’s accuracy and efficiency. 

II. Formulation of Optimal Approach 

To formulate the optimal approach, major aspects of the 
line simulation will be analyzed, viz., the formulation, the line 
characterization , the line model , and, the transient simulation 
method (see Fig. 1). Alternative approaches in each of the 
aspects will be compared to find the optimal combination 
that results in the maximum efficiency , accuracy and practical 
applicability. 

A. Formulation 

The formulation affects dimensions of the problem. One 
can distinguish between time-and-space formulations and time- 
only formulations. 

Time-and-space formulations (such as segmentation models 
[3], [4]) are based on the voltage and current distributions 
inside the line. These formulations are multidimensional and 
computationally extensive. 

Time ‘Only formulations deal exclusively with the voltages 
and currents at the line terminals. These formulations are one- 
dimensional and more efficient. Consequently, to achieve the 
maximum efficiency, the line simulation should be based on 
a time-only formulation. 

B. Line Characterization 

As can be observed from the system diagram [5] shown 
in Fig. 2, a line with terminations forms a feedback system. 
Therefore, one can distinguish between closed - and open-loop 
characterizations. 

Closed-loop characterizations (such as Z-, Y-, H - and S- 
parameter characterizations [6], [7]) include reflections from 
the terminations and lead to complicated oscillating transfer 
functions and transient characteristics (unit-step responses). 

The open-loop characterization (direct characterization in 
terms of the propagation functions) separates forward and 
backward waves and results in the simplest transfer functions 
and transient characteristics (see Fig. 3). The complexity of the 
transfer functions and transient characteristics is an important 
factor affecting accuracy and efficiency of the transient sim- 
ulation. Consequently, to attain the maximum efficiency and 
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Fie. 1. Aspects of the transmission line simulation. 



Fig. 2 . System model for transmission lines. W V f and W V b represent 
the forward and backward matrix propagation functions for voltage waves: 
T V i Tva an d IVi - IV2 stand for the near ’ ^ far ’ end matnx transmis- 
sion and reflection coefficients. 1 2 

accuracy, the line simulation should be based on the open- 
loop characterization. The complete set of expressions for the 
open-loop functions for uniform lines is given in Appendix A. 

C. Line Model 

Line models can be divided into two large groups: circuit 
and noncircuit 

Noncircuit models can not be directly integrated into a 
circuit simulator. As a result, these models cannot be effi- 
ciently applied the transient analysis of real circuits containing 

1 Throughout the paper, capital boldface, small boldface and normal italic 
symbols denote matrices, vectors and scalars, respectively. 

2 Since only muluconductor lines will be considered, the modifier "matrix” 
will be omitted in the future for brevity. 


thousands of nonlinear active devices. Examples of noncircuit 
models are the scattering-parameter model [6] and system 
model shown in Fig. 2. 

Circuit models can be directly incorporated into a circuit 
simulator and are of prime practical interest. They relate 
voltages and currents at the line terminals and do not depend 
on the terminations. Circuit models can, in turn, be subdivided 
into equivalent-circuit and device models. 

Equivalent-circuit models have a larger number of nodes 
than the line they represent. Examples of equivalent-circuit 
models are lumped and pseudo-lumped segmentation models 
[3], modal decomposition models for muiticonductor lines [8], 
[3], [9], and equivalent-circuit modeling of the propagation 
function and characteristic impedance based on Pade synthesis 
[10]. 

Device models have the same number of nodes as the line 
they represent. A well-known example of device models is 
the method of characteristics [1 1]— [15]. The circuit simulation 
time is cubically proportional to the number of nodes and to 
the number of voltage and current variables. Consequently, 
to achieve the maximum efficiency and practical applicability, 
the line simulation should be based on a device model that 
does not require the introduction of current variables. 


D. Transient Simulation Method 

The transient simulation method is the prime factor affecting 
the efficiency of the line simulation. The selection of the 
transient simulation methods is confined to the numerical 
Fourier and Laplace transformations \ numerical convolution 
an ^numerical integration. 
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Fig. 3. Examples of (a) open- and closed-loop transfer functions and (b) 
transient charactenstics. The scattering parameters correspond to the matched 
lossless relerence system. 

For the transformation and convolution methods , the com- 
putational complexity is higher than linear. These methods 
cannot directly handle nonlinear and time-varying systems, 
and, along with the discretization error, introduce time- and/or 
frequency-response truncation errors. In addition, the trans- 
formation methods cannot be directly used with recursive 
time-domain solvers employed by circuit simulators, that leads 
to relaxation techniques which additionally degrade the overall 
efficiency and accuracy [12]. Examples of the convolution- 
based techniques are the Spice model for lossy lines [9] and 
the scattering-parameter approach [6]. 

Numerical integration has linear computational complexity; 
it can directly handle nonlinear and time-varying systems, 
does not introduce truncation error, supports variable time- 
stepping and is compatible with recursive time-domain solvers. 
Numerical integration methods can be subdivided into direct 
and indirect. 

Direct numerical integration is based on approximations 
for integrals or derivatives and includes such conventional 
methods as linear multistep formulas, Euler, Euler-Cauchy and 
Runge-Kutta techniques. 


Indirect numerical integration [16] is based on the time- 
response invariant discrete synthesis, and has ideal accuracy, 
convergence and stability properties. Consequently, to achieve 
the maximum efficiency, accuracy and practical applicability, 
the line simulation should be based on indirect numerical 
integration. Indirect numerical integration covers as special 
cases techniques such as recursive convolution [13], [17], 
and approximation of the response of a linear network to an 
arbitrary piecewise linear input waveform used by some of the 
asymptotic waveform evaluation (AWE) methods [18]. 

To systems characterized with samples of their responses, 
numerical integration is applied via the difference approxima- 
tion method [19]. The method is based on the approximation 
of the system response with the corresponding response of a 
system for which numerical integration formulas are already 
available. The complexity of the difference approximation 
method is only that of the approximation itself. As soon 
as a system response has been approximated, the numerical 
integration formulas are readily available directly in terms of 
the approximation parameters. 

To attain the maximum efficiency and accuracy, the differ- 
ence approximation should be applied in the domain of the 
system characterization. For transmission lines it usually is 
the frequency domain. The time-domain approximation should 
be used only when time-domain responses are available. The 
complete set of analytical expressions for the fundamental 
time-domain open-loop responses of two-conductor uniform 
constant-parameter lines is given in Appendix B. It also 
includes a new simple and accurate asymptotic approximation 
for the responses of propagation functions. 

To improve accuracy, the delay should be separated from the 
propagation functions before the difference approximation is 
applied. The conventional frequency-domain method of delay 
separation for multiconductor lines is based on diagonalization 
with the frequency-dependent modal transformation matrices 
[14], [15], [17]. These matrices are nonminimum-phase func- 
tions of frequency with unstable time-domain responses, that 
limits the applicability of the modal transformation to special 
cases in which the matrices are constant [9], [20]. 

A novel matrix delay separation method [19] avoids the 
use of the frequency-dependent modal transformation matrices 
and is applicable to a general case of matrix transfer functions 
containing delay. For uniform lines, the formulas for the matrix 
delay separation from the propagation functions are included 
in Appendix A. 

£. Approximation Methods 

The choice of the approximation method for the difference 
approximation affects the overall efficiency, accuracy and reli- 
ability of the line simulation. Based on approximation cntena, 
approximation methods can be categorized into four major 
groups: minimum maximum error based methods , least square 
based methods , interpolation methods , and series expansion 
based methods (see Fig. 4). 

Mini-max methods provide the highest accuracy, but result 
in the most inefficient and unreliable algorithms (nonlinear 
optimization). 
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Fig. 4. Approximation methods. 


TABLE 1 

Comparison of Approximation Methods 


Approximation Method 

Maximum 
Relative Error 

Relative 

runtime 

Mini-max approximation 

0.0004% 

200 

Least squares approximauon 

0.002% 

50 

Interpolation 

0.003% 

1.0 

Economized rational approximation 

0.004% 

2.0 

Pad 6 synthesis at zero 

0.06% 

1.3 

Padd synthesis at infinity 

0.3% 

1.5 

J 


Least square methods provide high accuracy, but are still 
computationally extensive. The examples of least squares 
based methods for the time-domain difference approximation 
are Prony’s method [21] and pencil-of-funchon method [22], 

Series expansion based methods (such as Pad* synthesis 
used for the AWE [23]-[25]) are computationally efficient, 
but provide the poorest accuracy. 

Interpolation (point-fitting) [19] agrees exactly with the 
original function on a given set of samples. It provides high 
accuracy for simple functions, such as open-loop transmission- 
line responses, and is the most efficient among the approxima- 
tion methods. It also requires the minimal number of the orig- 
inal function samples, which is important when the samples 
are obtained from electromagnetic simulations. Consequently, 
to achieve the maximum efficiency, the line simulation should 
be based on interpolation. 

Table I presents the typical values of the maximum relative 
error in the full frequency range from zero to infinity and 
relative runtime for various approximation methods as applied 
to the third-order frequency-domain difference approximation 
of open-loop transmission-line functions. Economized rational 
approximation starts with the Pad* synthesis, which is not 
accurate away from the expansion point, and then perturbs it to 
reduce the leading coefficient of error in a given approximation 
interval [23], [24]. 

As one can observe, interpolation provides accuracy com- 
parable with that of the least square approximation, and 


is 200 times more efficient than mini-max approximation. 
Interpolation is also up to two orders of magnitude more 
accurate and 30-50% more efficient than Pad* synthesis. 
Note also that, because of the simplicity of the open-loop 
characterization, a very high accuracy is achieved with as low 
as third-order approximation and as few as seven samples of 
the original function used for the interpolation. 

F. Summary of Optimal Approach 

To summarize, the analysis of the problem showed that 
to achieve the maximum efficiency, accuracy and practical 
applicability, the line simulation should be based on 

• time-only formulation, 

• open-loop characterization; 

• device model that does not require the introduction of 
current variables; 

• indirect numerical integration; 

• frequency-domain difference approximation based on the 
interpolation and matrix delay separation. 

A method close to the optimal, but based on the time- 
domain difference approximation, was first proposed by Sem- 
iyen and Dabuleanu [17], and was further developed by 
Gruodis and Chang [15] to accommodate the frequency- 
domain approximation. The advantages of the approach were 
recognized only in recent years, and an increasing number 
of techniques close to optimal are published [13], [14], in- 
cluding techniques based on the AWE, recursive convolution 
and method of characteristics, and by researchers previously 
advocating transformation and direct convolution techniques. 
The applicability of the methods, however, had been limited 
by the lack of accurate, reliable and efficient frequency- 
domain approximation and delay separation techniques such as 
interpolation-based approximation methods and matrix delay 
separation [19], and by the lack of open-loop models for 
nonuniform lines. 

The authors’ implementation of the optimal method for uni- 
form lines is outlined in the next section. The implementation 
of the method for nonuniform lines is described in [26]. 

III. Implementation of Optimal 
Method for Uniform Lines 

A. Frequency-Domain Line Model for Transient Analysis 

The frequency-domain element characteristic (for the tran- 
sient analysis) which does not require the introduction of 
current variables and is suitable for the line modeling, is given 
by 

ri 1 M = YiMv 1 (w)-jiM 

\ l2(u>) = Yj(«)v 2 (w) - j2(w). 

The conventions for the terminal voltages and currents are 
shown in Fig. 5. The expressions relating the matrix admit- 
tances Yi and Y 2 and vector current sources ji and j 2 to the 
transmission line characteristics are derived directly from the 
continuity conditions for the voltages and currents at the line 
terminals. To separate forward and backward waves and open 
the feedback loop, the current source ji must depend only on 
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Fig. 5. Conventions for the voltages and currents at the line terminals. 


the backward wave, and j 2 only on the forward wave. This 
condition uniquely defines Yi,Y 2 and ji,j 2 as follows: 

YiM = Y 2 (w) = Y c (u/) (2) 

and 

f jiM = 2ibi(^) 

\ ja(^) = 2ifj(u/). 


Once the approximation has been performed, indirect nu- 
merical integration formulas (discrete-time difference equa- 
tions) are readily given in terms of the approximation param- 
eters. For the step invariance the formulas are 

f Af 

i y(tn) = HooXitn) + ^ r m («„) 

m = 1 

t z m (t n ) = a m (l - e—‘- T “ )x(t„_i) + Zm(tn-l)- 

( 6 ) 

where x, y and z m stand for the excitation, response and state 
variables, respectively, and T n — t n - t n - 1 is the time step 
at the nth transient iteration. 

For the ramp invariance 

r M 

< y{tn) ~ Ho %{tn) ~ ^ (fn ) 

m= 1 

l = d m (T n )(x(t n ) - x(t„-i)) + e-“— r "c m (tn-i). 

(7) 


where Y c stands for the characteristic admittance, and the 
forward and backward current waves, in, in and ibi,ib 2 are 
related as follows: 

f ibiM = Wib(u;)[ib2(^) = i2(^) + ir2(^)] 

\ in(tu) = Wif(cj)[in(^) = ii(t^) + ibi (<*>)]• 

For uniform lines, the propagation functions for the forward 
and backward current waves are equal, Wjf = Wib- The 
propagation function and characteristic admittance can be 
computed from the insertion loss data [15], scattering param- 
eters [27], or distributed RLGC parameters (see Appendix 
A). 

As can be observed, for uniform lines, the open-loop device 
model ( 1 )— (4) is equivalent to the generalized method of 
characteristics [12], [13], [15]. However, for nonuniform lines, 
the generalized method of characteristics no longer separates 
forward and backward waves and loses physical meaning [26]. 


B . Difference Approximation 

To perform the transient analysis, indirect numerical in- 
tegration [16] is applied to the propagation functions and 
characteristic admittances in the frequency-domain line model 
(IM4) by using the difference approximation method [19]. 

For the difference approximation in the parallel canonic 
form, samples of the frequency-domain transfer function are 
approximated with the rational polynomial function 


A/ 


ff(jw) = H ao +Y,J 


m= 0 


+ jw/uv 


(5) 


or samples of the time-domain unit-step response are approx- 
imated with the exponential series 


M 

h(t) = Ho ~ 

m= 1 


where H 0 and H 0 c denote the initial and final values of the 
approximating transfer function H(juj). 


where 

^ 1 - e - ~' rmTn 

dm{Tn) — — ■ 

^cm-L ri 

An alternative form of the ramp-invariant indirect numerical 
integration formula has the coefficients of the present-time 
sample of the excitation lumped together 

r / m \ a/ 

y{t n ) = ( Ho - ^2, ] X (fn) - ^2 

\ m= 1 / m= 1 (8) 

z m (t n ) = (dm(T n .i)e— - r - - d m (T n ))x(tn- 1 ) 
z m (t n - 1 ). 

This form is especially suitable for discretization of charac- 
teristic admittance, because, for admittances, the present- and 
past-time terms of the numerical integration formulas have 
different physical meanings. 

Before the approximation, the delay is separated from the 
matrix propagation function using the matrix delay separation 
formulas from Appendix A, and is modeled separately using a 
low-order spline of the simulated time points. The difference 
approximation is applied to each element of the delay less 
propagation function and characteristic admittance matrices. 
For the characteristic admittances in (1), the excitations are the 
terminal voltages and the responses are the terminal currents. 
For the propagation functions in (4), the excitations and 
responses are the current waves. 

Since the open-loop functions are aperiodic, they have to 
be approximated with only real poles, — u; rm . In addition, the 
poles have to negative to be assure stability. 

To represent the original functions accurately with the 
minimum number of samples, the variation of the original 
function from sample to sample should be about the same 
The following empirical formula for the sampling frequencies 
was found to provide good results 


Uk 




k = 0.1. - .A". 


The end of the approximation interval, w* , should be cho- 
sen so that the original function would closely approach its 
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final value. This assures that the resulting indirect numerical 
integration formulas will be accurate in the full frequency and 
time ranges from zero to infinity. 


C. Interpolation-Based Complex Rational Approximation 
Method for Frequency- Domain Difference Approximation 

The method fits samples of a complex transfer function 
H{ui) with the rational polynomial function (5) at the set of 
arbitrary spaced frequencies (0. ui . • • • . u>k }. The method 

proceeds in three steps. 

First, the real part of the original function is fit with the real 
part of the complex rational polynomial function, which is a 
real rational polynomial function of squared frequency [19] 


Re{H(juj )) 


c 0 + cioj 2 4- C2(or) 2 + • • • + CAr(u>~) A/ 

1 + + 0 2(w 2 ) 2 + • • • + iht(u 2 ) sl 


The following linear system of equations (10) (see bottom of 
the page) results from matching the real part of the original 
function with (9) at the set of frequencies and premultiply- 
ing both sides of each equation with the denominator For 
interpolation, K = 2M and solving (10) produces a rational 
polynomial function which coincides with the real part of the 
original function at all of the sampling points. For a set of 
samples larger than 2M + 1. the least square solution of (10) 
can be obtained. However, it minimizes the approximation 
error premultiplied by the denominator, which can lead to 
inaccurate approximation. Better results are achieved with the 
method of averages [28], which partitions the larger number 
of equations into 2A/+ 1 subsets in the order of the increasing 
of w. The equations within each subset are added up, which 
makes the system consistent. The method is effective in 
averaging out the noise in measured data. 

After the real part has been approximated, the denominator 
of (9) is factored yielding the squared poles, -a J 2 m . Conse- 
quently, no unstable right-half-plane poles can be produced. 
However, there still can be spurious complex conjugate and 
purely imaginary poles, which are removed. The remaining 
real negative poles are used to formulate the equations for the 
partial expansion coefficients, a m , of (5). As a result, the order 
M of (5) is less or equal to that of (9). 

Matching the real and imaginary parts of the onginal 
transfer function ff(u) with the corresponding parts of (5) 
at the set of frequencies {0,uti,u> 2 , ■ ■■ ,^k} leads t0 the 


following linear system of equations 
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( 11 ) 


Um(H(u K ))l 


For interpolation, M = 2 K. and both real and imaginary parts 
of the original transfer function are matched exactly at all of 
the K frequency points and dc. For an arbitrary larger number 
of points, the least square solution of (11) is obtained from 

A t Ax = A r b. 

The total computational complexity of the approximation 
method is that of two real linear solutions and one polynomial 
factoring. The orders of the polynomial and linear systems 
depend only on the order of the approximation and not on the 
number of the original function samples. Since no iterative or 
relaxation techniques are involved, the method is free from 
convergence problems. The method can be extended to match 
exactly the initial and final values of the original function and 
to perform a complex-pole approximation [19]. 

Fig. 6 shows an example of the fourth-order approximation 
of an open-loop transmission-line transfer function. As can be 
observed, although only nine samples of the original function 
were used, the approximation exhibits an excellent match in 
the full frequency range. In general, due to their simplicity, the 
open-loop functions can be accurately fit with the 3rd-9th- 
order approximation. 
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Fig 6. An example of the 4th-order complex rational approximation. The 
original function is shown by the thin continuous curves and the approximating 
function by the thick dashed curves 


D. Companion Model 

By applying the difference approximation to the propagation 
function and characteristic admittance, the frequency-domain 



element characteristic (1) is transformed into the following 
discrete-time element characteristic, or companion model 

f il(tn) = Yi(* n )Vl(t n ) ” JlUn) (12) 

lh(U = Y 2 (^)v 2 (U-j2(<»)- 

The circuit-diagram interpretation of the companion model is 
shown in Fig. 7. 

The admittances Yi and Y 2 represent present-time co- 
efficients in the indirect numerical integration formulas for 
the admittances Yi and Y 2 . The current sources ji and j 2 
combine the currents jYi and jYr which correspond to the 
remaining parts of the numerical integration formulas for the 
admittances, and ji and j 2 , which are given by the discretized 
(3) and (4) 

/ jl(*n) = -jYi(*n) + JlUn) (13) 

1 j 2 (*n) = “JY 3 (^n ) + JsUn)- 

Equations (3) and (4) do not contribute to the admittance part 
of the companion model because the propagation functions 
contain a delay. 

The Modified Nodal Approach (MNA) stamp corresponding 
to the companion model (12) is (see (14) at the top of the 

page). . 

In the circuit simulator during the transient analysis, tne 
lines are represented by the tables of numbers (14), which 
are recursively updated at each time iteration using numerical 
integration. The left-hand side of the stamp (14) has to be 
updated only when the value of the time step changes. If the 
step-in variant indirect numerical integration formulas (6) are 
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used, the LHS of the stamp becomes independent of the time 
step, and only the right-hand side vector has to be updated. 

Since the terminal currents are not introduced as variables, 
the values of ii and i 2 in (4) are computed from (12). 

E. Line Model for AC and DC Analyses 

For ac and dc analyses, the complexity of the transfer 
functions is not important, and the element characteristic that 
does not require the introduction of current variables and is 
suitable for the ac/dc modeling of transmission lines is the 
V-parameter characteristic 

f iiM = Yu(u/)vx(u/) + YtaMvaM (15) 

\ hM = YaxMvxM + Y22 (w) v 2(w). 

The Y'-parameters are related to the open-loop functions as 
follows: 

Yn(w) =Ya2(u/) 

= Y C M + 2 [I - W* (w)] _l Wj (w)Y c (u»). 
Yia(w) »= YaxM = -2(1 - W?(ai)] _1 Wi(w)Y t («), 

where I is the identity matrix. The expressions were derived 
by eliminating j i and j 2 from (1)— (4), and transforming them 
to the form of (15). 

The dc model is merely the ac model at zero frequency. For 
the limiting case of lines with zero distributed conductance, 
G = 0, the dc values of the Y -parameters are 

Yn(0) = Y 2 a(0) = -Yia(O) = -Y 2 i(0) = y R_1 - 

The MNA stamp corresponding to (15) is (see (16) at the 
top of the page). 

F. Initial Conditions for Transient Analysis 

The dc model is used to perform the operating-point (op) 
analysis before the transient simulation. The op solution is 
then used as the initial conditions for the transient analysis. 
The initial conditions for the indirect numerical integration are 
the dc values of the state variables, which are related to the dc 


value of the excitation, xq. as follows: z m (t o) — for the 

step-invariant case (6). z m (t») = 0 for the ramp-invanant case 

(7) , and r m (< 0 ) = -d m (Ti)x 0 for the ramp-invanant case 

(8) . The dc values of in and i b2 . which serve as excitations 
for the propagation functions in (4), have to be expressed in 
terms of the terminal voltages obtained from the op analysis. 
Resolving (1M4) leads to 

f in(0) = (I — Wj (0)] -1 [Y c (0) Vl (0) 
-Wi(0)Y c (0)v 2 (0)] 
i i b2 (0) = [I — Wj (1)] -1 [Y c (0)v 2 (0) 
-W,(0)Y e (0)vx(0)]. 

For the limiting case of G = 0. the expressions become 

irx(0) = — i b2 (0) = ^R- _1 [ v i(°) - v 2 ( 0 )]- 

G. Optimal Line Simulation Algorithm 

For an MNA-based simulator, the optimal line simulation 
algorithm is as follows: 

1) Before the transient analysis: 

a) Perform op analysis of the circuit to find the initial 
conditions for the transient analysis. Use the ac/dc 
model (15M16). 

b) For each line in the circuit, perform the difference 
approximation of each element of the propagation 
function and characteristic admittance matrices. 

2) At each time iteration: Recursively update the line 
stamps using the indirect numerical integration formulas 
obtained at step 1(b) and companion model (12)— (14). 

Since the method introduces neither additional nodes nor 
current variables, the optimal line modeling does not increase 
at all the circuit solution time. The only additional time is 
required to perform a low-order interpolation once in the 
beginning of the simulation, and for a low-order numerical in- 
tegration. As shown in the next section, this time is negligibly 
small compared to the circuit solution time. 
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TABLE II 
Relative Runtimes 





Relative Runtime* 


Number 

of 

Nodes 


Lines 

Lines Modeled with the Optimal Method 

Circuit Description 

Modeled 

with 

Total 

Indirect 

Difference Approximation* 1 


Lumped 

Resistors 6 

Circuit 

Simulation 0 

Numerical 

Integration** 

Frequency 

Domain 

Time 

Domain 0 

10 

2 two-conductor lines. 1 four -conductor line 
(lines and excitation sources only) 

0.937 

LOO 

0.0927 

0.0793 

0.0550 

10 

2 two-conductor lines, 12 MOSFETs 

1.00 

1.00 

9.66 10' 3 

9.47 10" J 

6.57 10’ 3 

100 

20 two-conductor lines, 10 four-conductor 
lines (lines and excitation sources only) 

0.999 

1.00 

2.06 10' 3 

1.75 10"* 

1.21-10'* 

1000 

200 two-conductor lines, 100 four-conductor 
lines (lines and excitation sources only) 

1.00 

LOO 

1.64-10"* 

1.39 10"’ 

9.64 10“ 


• For 1000 time points. 

b One resistor per signal conductor. 

c Does not include the difference approximation time. 

d Seventh-order 

e Time-domain difference approximation was performed by the relaxation interpolation method [191. The runtime includes 
automatic determination of the approximation interval and interpolation points. 


H. Numerical Results 

The optimal method has been adopted in several industrial 
and commercial circuit simulators, and, in numerous real- 
life simulation exercises, proved to be reliable, accurate and 
efficient. Table II presents relative runtime data for circuits of 
various types and sizes. As can be observed, even for the worst 
case of a small circuit consisting only of lines, the optimal 
model is virtually as efficient as the simple replacement of 
interconnects with lumped resistors. The resistive model was 
chosen for the comparison because it represents the limiting 
case in the simplicity and computational efficiency of the 
interconnect modeling. 

Fig. 8 shows verification of the optimal model accuracy with 
the Spice3e2 lossy multiconductor line model [9]. A simple 
network was chosen as an example to reduce the influence of 
factors other than the line model on the simulation accuracy. 
A variable, third-to-fifth-order frequency-domain difference 
approximation was applied. As can be observed, the compared 
waveforms are indistinguishable. In fact, the accuracy of 
the optimal method depends exclusively on the accuracy of 
the difference approximation which is very high (see Table 
I). The runtime for the optimal model was three orders 
of magnitude shorter than that for the Spice model, which 
is, in turn, an order of magnitude faster than segmentation 
models. 


IV. Conclusion 

From a novel analysis, based on identification of the most 
significant aspects of the problem and comparison of alter- 
native approaches in each of the aspects, it was shown that 
to achieve the maximum efficiency, accuracy and practical 
applicability, the line simulation should be based on: a time- 
only formulation, open-loop characterization, device model 
that does not require the introduction of current variables, in- 
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Fig. 8. (a) The network configuration and (b) comparison of the transient 

waveforms generated using the optimal line model installed in an MNA-based 
circuit simulator (thick broken curves) and Spice3e2 (thin continuous curves). 
/?, = /?.-, = 50 Si./?* = /?<, = 1 k Si./?:, = Ra = MU: 
self-inductance L s = 416 nR/m. self-capacitance C\ = 04 pF/m. mu- 
tual inductance L m = 125 nH/m, mutual capacitance Cm = 22 pF/m. 
R — 15 Q/m.G = ()./ = 0.G77 m (all signal conductors are the same). 

direct numerical integration, and frequency-domain difference 
approximation based on the interpolation and matrix delay 
separation. 
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TABLE ill 


Name of Function 

Function for Current Waves** 

Funcuon for Voltage Waves 4 

In Standard Basis 

In Modal Basis 

In Standard Basis 

In Modal Basis 

Admittance and 
impedance per unit 
length 

Y(oi) = G(o>) + ;<i jC((D) 
= M,(co) Y„(co) M» (ai) 

Y.«d> 

= Mf'(o» Y((i)>M v (a)) 

Z(io) = R(u)) + ;(i)L(tu) 
= M v (o» Z„(w) M7'(o» 

Z.(o» 

= My (tll)Z((l>) M,((0) 

Modal 

transformation 

matrices 

M,(co) 

= Eigenvectors^ (to) Z(<o)), 
M, = M,(~) 

= Eigenvectors^ C(<*>) L(°°)) 

Mi'(a>) = Mj(ai), 
M,- 1 = Mj 

M v (d)) 

= Eigenvector^ Z(<i>) Y(o>)), 
M V = M V («) 

= Eigenvectors(L(«*) C(«»)) 

My (01) = M*(ui), 
My‘ = Mf 

Characteristic 
admittance and 
impedance 

Y c ((o) = K,«i))Z-'(o)) 
= Kj'(w) Y(u» = Z;'(ti)) 
= M,(ci>) Y„(o»Mv(co) 

1 

Y_(o>) = (Z''(ai) Y.(w))» 
= K^((i>)Y.(o)=Z^(u)) 
= M;'«o) Y,(co) M v (oj) 

Z,(oi) = K v (oj) Y '«d) 
= K,(a») Z(u)= Y;‘((0) 
= M v (a>) Z — (oi) Mj'tco) 

Z < .(oi) = (Y; , (o»Z.(o»))j 
= K;'.((1))Z.((11) = Y^(<J1) 
= My (oi) Z,(oi) M,(oi) 

Propagauon 

constant 

K,(ui) = (Y(o» Z(u))i 
= Y(cd) Y;'(o)) 

= M 1 (co)K t .(a))Mf'((o) 

K,.(a)) = (Y.(a))Z.(m))i 
= Y.(oj) Y^(oo) 

= Eigen vnlue^K, (to)) 

= K Vb ((0) 

K v (oo) = (Z(a>) Y(cd))i 
= Z(o))Z t "‘(ti)) 

= M v (co) K B (to) M;‘(<d) 

K v .(®) = (Z il (w) Y.(o)))> 
= Z B ((0) Z^(co) 

- Eigenvalues^ (tu)) 
k v «o)M v «i»= 
= K ta (o» 

Propagation delay 6 

T, =(C(») L(-))'» l 
= K,(~M = M, T ta Mf' 

T ta = K,_(~) / 

= Eigcnvaluc^T, ) 

*m;‘ t, m, =t v . 

T v =(L(-)C(-))i/ 

= K,(-)/ = M, T v _ My 1 

Ty.=Ky.(-)( 

= Eigenvalue^ T v ) 
= m; 1 T v M, = T_ 

Propagation 

function 6 

= Y,(ai)W v (o» Z,(o» 
= W,(a>) M, M; 1 

W to (u) = 

= W,(ui) M,(a>) 

= W v .(to) 

W v (u>) = e- ,t '’ ,, ■ ,, 

= M v (o»)W^(to) My (oi) 
= Z,(o» W,(o» Y,(ai) 

= W v (a>)e' J-T ' 

= Wy(Ol) My e"'“ T — My 

W v .(oi) = 

= m;‘(oi) w v ((o)My(oi) 
= W„(oi) 

Delayless 
propagauon 
function 6 

W I (to) = W,(to) t"*' 


W v (<o) = W v (o))e'“ T * 


1 

Transmission 

coefficient 0 

T,(co) = (I + Y,(w) Z,(a>)) ' 
= Y,«d)T v ((o)Z,(m) 


Ty(ai) = (I + Z,((o) Y.(ai))" 1 
= Z,«o)T,(o>)Y,<a» 


Reflection 

coefficient 0 

r,(a>) = (l + Y,(<o) Z,(o>)) ' 
(i - Y,(ai) Z,(oi)) 

= Y e (ai)r v (oi)Z c (a>) 


r v (oi) = (Z,(oi) Y,(a>)+ 1) 1 
<Z,(a» Y.(o»-I) 

= Z,(o» r,(oi) Y c (o» 



* Voltage and current functions are related via the following duality replacement rules: 

i 

b Boldface (.)* and e u denote matrix squrc root and matrix exponential, respectively. 
c For a Thevenin's termination Z^to). 


The practical implementation of the optimal method for 
multiconductor lossy frequency-dependent lines characterized 
by discrete samples of their responses was outlined, including 
extraction of initial conditions from op analysis and the line 
model for ac/dc analysis. The complete set of expressions for 
the open-loop transmission-line functions was given, including 
new formulas for the matrix delay separation from the propa- 
gation functions, which avoid the use of frequency-dependent 


modal transformation matrices. The complete set of analyt- 
ical expressions for the fundamental open-loop time-domain 
responses of two-conductor lines was presented, including a 
new simple and accurate asymptotic approximation for the 
responses of propagation function. The novel interpolation- 
based complex rational approximation method was introduced, 
and ramp- and step-invariant indirect numerical formulas were 
given in terms of the approximation parameters. 
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TABLE IV 


Name of 
Function 

Propagation 
fun ebon, 

W( o» a 


Character* sbe 

admittance. 

K r (u» c 


Transmission 
coefficient for 
voltage waves, 
T v «o) c - d 


Arbitrary 

G 


MO 


G = 0 


Arbitrary 
C 


G = 0 


Arbitrary 
G 


Reflection coe- 
fficients r(( 0 )» 


G = 0 


Transient Characterisbc, h(t) 




MO - ( 1 - - 7 = “(' - V b 

y2nb(i + d) 


hr(l ) = T 0 + § J<- f ’/o(*T) ) 


Impulse Characteristic, g(») 


M») = *-*6(0+ ae-" 




MO -^8(0+ (2jtt(f+(f)) l 


u(f -t) b 


grfl) = T 0 |6(/) + + <* /i(W | 


h,U)*Y t e*‘l t $t) 


w,)= ^ {-*[£♦ O-SH 

+ R, e-*‘Ubt) + OP -c)ir" JeO-o'/oO*) 

L 0 

i < 

+ ST /„(**) rfx d* 


+ * 


m,)= 

i 

c*'l<t(bi) + OP ~ c) e* W*t) 


gr (/) = T. {5(0 + P ^<"[/,(pt)-/o(P')] } 


04- ^| 2 «( c -f) 

+ /?, [e" p, [(P -c)/o(W + bl\{bt)\ 

sy-j- 


+ ^c(c-2P) + 




+ /?, [(P - c)/ 0 (W + b / i(W] 

r - ‘ 

+ c(c - 2P) hOn) dt 


M') = I -2A r (t) 


*K0= 8(0- 2grO) 


Note: «(,) is the unit-step function, a= ^GZo+ « = 2 | GZ ® “ Zo! ^ ~ 2 ( C + L )’ fr_ 2 |c L |’ 

_ K d ** t= >/aT/. T„ = a/J". and Z 0 = -j§". For C = 0: a = fl= ^- o and P = t= 2t- 

c_ L - Zo 2n«l 

» In the case of two-conductor lines, this functions are the same for both voltage and current waves. 

b Approx, mation based on the asymptodc expans, on of the modified Bessel funcuon /,(*>. h is within 1% accurate ,n the full tune 

^Expressions for the cooesponding dual (voltage/current) functions can be obtained via the duality replacement rules (see 
Appendix A). 

d For a Thevenin’s terminauon R,. 


The optimal method is compatible with recursive time- 
domain solvers employed by circuit simulators and supports 
variable time-stepping. The method has been adopted in sev- 
eral industrial and commercial circuit simulators, and, in 
numerous real-life simulation exercises, proved to be reliable 
and accurate. It was shown on an extensive set of runtime data 
that, based on the optimal approach, accurate line modeling in 
a circuit simulator is as efficient as a simple replacement of 
interconnects with lumped resistors. 


APPENDIX A 

Open-Loop Transmission-Line Transfer Functions 
See Table 111. 


Appendix B 

Time-Domain Open-Loop Responses of 
Constant-Parameter Two-Conductor Lines 

See Table IV. 
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Multilayered Dielectric Medium Using 
Closed-Form Spatial Green’s Functions 

Kyung Suk Oh, Student Member , IEEE, Dmitri Kuznetsov, and Jose E. Schutt-Aine, Member . IEEE 


Abstract — An efficient method to compute the 2-D and 3- 
D capacitance matrices of multiconductor interconnects in a 
multilayered dielectric medium is presented. The method is based 
on an integral equation approach and assumes the quasi-static 
condition. It is applicable to conductors of arbitrary polygonal 
shape embedded in a multilayered dielectric medium with possi- 
ble ground planes on the top or bottom of the dielectric layers. 
The computation time required to evaluate the space-domain 
Green’s function for the multilayered medium, which involves 
an infinite summation, has been greatly reduced by obtaining 
a closed-form expression, which is derived by approximating 
the Green’s function using a finite number of images in the 
spectra] domain. Then the corresponding space-domain Green’s 
functions are obtained using the proper closed-form integrations. 
In both 2-D and 3-D cases, the unknown surface charge density 
is represented by pulse basis functions, and the delta testing 
function (point matching) is used to solve the integral equation. 
The elements of the resulting matrix are computed using the 
closed-form formulation, avoiding any numerical integration. 
The presented method is compared with other published results 
and showed good agreement Finally, the equivalent microstrip 
crossover capacitance is computed to illustrate the use of a 
combination of 2-D and 3-D Green’s functions. 

I. Introduction 

I N recent years, the characterization of microstrip disconti- 
nuities in a multilayered dielectric medium by equivalent 
circuits has gained special interest due to the modem develop- 
ment of VLSI technology. For an inhomogeneous structure 
the modes are hybrid, and the full-wave analysis must be 
considered. However, the quasi-static approximation is suffi- 
ciently correct when the transverse components of the electric 
and magnetic fields are predominant over the longitudinal 
components; in other words, the transverse dimensions of 
microstrip lines are much smaller than the wavelength. Based 
on this quasi-static approximation, we present an efficient 
method to compute the 2-D and 3-D capacitance matrices 
of multiconductor interconnects in a multilayered dielectric 
medium. 

Under the quasi-TEM approximation, the capacitance cal- 
culation follows from the solution of Laplace’s equation 
with appropriate boundary conditions. Various methods have 
been employed to obtain the solution in two-dimensional 
space [l]-[9]. TWo commonly used techniques for both 2-D 

Manuscript received April 5, 1993; revised October 11, 1993. 

The authors are with the Electromagnetic Communication Laboratory, 
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and 3-D capacitance calculations in multilayered structures 
are the integral equation method [10], [11] and the domain 
methods, such as the finite difference and finite element 
methods [12], [13]. In the domain methods, the unknown 
potential distribution is solved to compute the capacitance 
over an entire domain by either directly approximating the 
differential equation with the finite difference equation (FD) or 
using the equivalent variational expression in conjunction with 
the method of subareas (FEM). The major disadvantage of the 
domain methods is that the unknown potential distribution to 
be sought is over the entire geometry considered, including the 
dielectric regions; hence, it may be computationally inefficient 
for the open geometry case even with the use of absorbing 
boundary conditions to truncate the open geometry. On the 
other hand, the integral equation approach first obtains the 
Green’s function for a layered medium using image theory, 
which consists of rather slowly converging infinite series, and 
solves for the charge density on the conductor surfaces using 
this Green’s function as its kernel. As noted in [2] for N layers, 
the expression for the Green’s function would consist of N - 1 
infinite series. Alternatively, the free-space Green’s function is 
used in [2], [3] to avoid infinite series, but additional unknown 
charges on the dielectric interface and ground planes, on top 
of the unknown charges on the conductor surface, must be 
included. Hence, the dimension of the resulting moment matrix 
is substantially increased. 

Yet another approach to avoid an infinite summation is to 
solve the integral equation in the spectral domain (SDA), 
where the Green’s function is in a closed-form expression; 
however, this approach can not be applied to general problems, 
e.g., conductor with a finite thickness. In this paper, the 
Green’s function for the layered medium is approximated 
in the spectral domain using exponential functions, which 
is equivalent to a finite number of weighted real images 
in the space domain. Although complex-valued exponentials, 
which are often used in a nonquasi-TEM analysis [14], can 
also be employed to reduce the number of weighted images 
[15], real-valued exponentials are sufficient to approximate 
for quasi-TEM applications, and it further avoids the use 
of expensive complex operations. Since the spectral-domain 
representations of the Green’s function for 2-D and 3-D cases 
are identical, the approximation is only performed once for 
both cases, and then the equivalent weighted images in the 
spectral domain are directly used to evaluate the Green’s 
functions in the space domain. 


00 1 8-9480/94 $04. 00 © 1994 IEEE 
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Conductor 
in Liyer 1 


Conductor 
in Layer 2 


(b) 


Fig. 1. (a) Cross-sectional and (b) planar views of possible configurations of 

multiconductors in a multilayered medium. The two figures are not related. 


II. General Statement of the Problem 

The general geometries of systems of multiconductors em- 
bedded in a multilayered medium are illustrated in Fig. 1; 
(a) shows a possible cross-sectional view while (b) shows a 
possible planar view. An arbitrary number N* of nonmagnetic 
dielectric layers are backed by two optional ground planes 
with possible top or bottom locations, and within these layers 
an arbitrary number N c of perfect conductors are placed 
throughout the layers with arbitrary orientations and possible 
discontinuities. The geometries of the dielectric layers are 
assumed to be uniform in the x- and 2 -directions, and the 
cross sections and planar geometries of the conductors can be 
arbitrary as long as their boundaries can be described with a 
piecewise linear function. 

The integral equation relating the electrostatic potential 
V{r) to the charge density p(r) is 

V (r) = / G(r,r')p(r')dr f ( 1 ) 

J Q 

where G(r. r') is the Green’s function for the multilayered 
medium, and Q denotes the surfaces or cross-sectional bound- 
aries of conductors for 2-D or 3-D problems, respectively. The 
capacitance can be computed by solving this integral equation 
for the charge density p(r') with various settings of voltages on 
the conductors. We will first concentrate on the determination 
of the Green’s function G(r,r'). 

III. Approximation of the 
Spectral-Domain Green s Function 

Consider a unit point charge located at the mth layer at 
(^o,yo 5 ^o) (Fig- 2). By definition, the 3-D Green's function 



Fig. 2. the geometric configuration used for determining the Green's func- 
tion. 


satisfies Poisson’s differential equation: 


V 2 G 3D (i,y ,2 | x 0 . j/o, 2o) = -S(x - x 0 )6(y - y 0 )6(z - zq) 

( 2 ) 

with the appropriate boundary conditions at the possible 
ground planes and dielectric interfaces. Noting that the 
dielectric medium is uniform in two directions, we can 
represent the Green’s function and the point source in the 
spectral domain in terms of its transforms in the x- and 2 - 
directions. The space-domain and the spectral domain Green’s 
functions are then related by 


G 3D (x,y,z | x 0 .yo, 2 o) 


1 r-r 00 f+00 

= / / dadPe-^-^-^^- 20 ') 

(27T j J — oc J— 00 

x G 3 D (a. y, 0 I x 0 , 2 /o- 2 o) 


G 3D (a.y,0 | x 0 ,y 0 ,z 0 ) 


(3a) 


/ + OO r+oc 

/ dxdze JQ(x ~ Xo)+jS(l ~ 2o) 
-00 J — oc 


x G 6V (x.y,z | x 0 . yo,z 0 ) 
(3b) 


where G 3£) (a.j/./3 | xo^yo-^o) is the 3-D spectral-domain 
Green’s function and a and 0 are the transform constants 
associated with the x and 2 directions, respectively. Then, (2) 
can be written as 

(sj/2 ” q2 - /3 2 ^G 3D (a, y, 0 I x 0 ,yo,2o) = - Vo)- 

(4) 

The general solution of the above equation is given by 


G 3£> (7 ; y|r 0 ) = 


* — rv 


+ B„ 


2£m7 


7 = \A * 2 + 

(5) 


where the first subscript m denotes the layer where the 
source is located (source point), while the second subscript 
n denotes the layer where the Green’s function is evaluated 
(observation point). The same expression can be obtained by 
Fourier transforming G 2D {x,y \ x 0 l yo) in the x-direction 
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with the transform variable 7. Furthermore, since the unknown 
coefficients A m , n and B m , n are to be determined using the 
boundary conditions in the y-direction only, it is easily seen 
that the expressions for G 2D and G SD must be identical 
under our Fourier transformations. Thus, in what follows, 
we will omit the superscripts for the spectral-domain Green’s 
functions. Applying the boundary conditions at the dielectric 
interface and ground planes, ( 5 ) can be written as 

6(7,1/ 1 ro) = (e- 7 ° v + f„, n+1 e- 2 ^" + ™) y > yo 

(6a) 

6 ( 7 , 2/ 1 r 0 )= ^ L (e +7y +f n , n _ 1 c +2 ^"-^) y< y 0 

(6b) 


where 


n— 1 


At.n = n 5 ^ + l 


j=r 


Aj n,n A m , m 1 1 — 1 

j=n + l 

c+ _ Thill 


Tjj-i 


3 ' J ~ l i-r 


( 7 ) 


( 8 ) 


f _ T jJ+1 +r j+l . j+ 2e 2 ^- d ^A 

1,J+1 i + r JJ+1 f ;+1 , J+ 2e 2 ^-^‘) 

1,3 1 i + r i , J _ 1 f > _ 1 , J _ 2 e 2 ^-*-^-‘) 


r 


*0 


g» ~ £ J 

€ t + £j 


T,, = 


2cj 

£i + £j ' 


( 9 a) 

( 9 b) 


Here. f jj+i is the generalized reflection coefficient, which 
is the ratio of the amplitudes of voltages at y = d n due to the 
image charges located above and below the jth layer. fjj+i 
takes the value of 0 or - 1 if the jth layer is a half space, or the 
(j + l)th layer is a ground plane, respectively. The unknowns 
to be determined are 4 + m and A~ <m . Using the facts that, 
at y = yo, (6a) and (6b) must be equal and that the normal 
component of the displacement field must be discontinuous by 


the charge density at y — y 0 , we can obtain expressions for 
4 m.m and 4 -. m as follows: 

A+. m = M m [e + ™ + (10a) 

A", m = A/ m [e~™ + f m . m+1 e- 2 ^~ + ™] (10b) 


where 

M m = fm.m 


if 


m, 77 i + l 


g + 2-y(d m _ i —dm ) 


(ID 


Now to approximate the Green’s function with respect to 7 
but independently of y and yo* we rearrange these expressions 
by factoring out all y and yo dependencies as follows in (12a) 
and (12b), shown at the bottom of the page, where 


n — 1 


jcui = M m f n , n+1 n 5 

■+ 

J.J + 1 

j=m 

Tl — 1 

^m.n, 2 = M m T n,n + l^ m,m- 

-n*. 



71 — 1 


*£„.3 = Aim I! 5 ^ + l 


>=m 


n— 1 


^m,n,4 = Mm f m,m — 1 


j-m 


m 


^m,n, 1 = M m T m,m + l 


j = n-f 1 

m 


*m.». 2= Mm n 


j = n + l 

771 

■^m,n,3 = ^mFm.m + l Tn.n- 

n 


j = n + l 

m 


^m,n,4 = ^mf n,n- 1 J"J 

S jj- 1* 


j=n+l 


( 13 a) 


( 13 b) 


The determination of the closed-form space-domain Green’s 
function can now be preceded by the approximation of the 
coefficient functions K± n i of the exponential terms. 

A physically intuitive approach to approximate the potential 
due to a charge in the layered medium is the use of a finite 
number of the weighted image charges in the homogenous 
medium, which is equivalent to approximating the coefficient 
functions K± n i with exponential functions. The equivalence 


6 ( 7 , 2/ 1 r °) = 2~ (*m,n, 1 e 7(v+V0-2d " ) + K, n ^ {+V ~ V0+ 2 ^- M) + tf+n, 3 e 7(_V+Vo) 
+ ^m,n, 4 e 7( ' V_yo+2dm -‘ ) ) 2/ >2/0 

6 ( 7 ,: V I ro) = 2^(^m,n,ie 7(v+v °- 2,1 ”' ) + 


(12a) 


+ tf-, n ,4e 7( - v - V0+2d "- l) ) y < 


yo 


(12b) 
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Fig. 3. The geometry used to test the approximated Green’s function. 

between the weighted image charges and exponential functions 
will be shown below in (15) and (16). In EM analyses, 
the complex-valued exponentials are often used for pole- 
zero modeling of signals, such as an electromagnetic-scatterer 
response. The least-square formulation of this exponential 
approximation results in nonlinear equations and can only 
be solved by iterative methods, such as gradient descent 
procedures or the Newton method. Due to the computa- 
tional inefficiency of these algorithms, some other suboptimal 
noniterative techniques are also proposed: the least-squares 
Prony method and the generalized pencil-of-function (GPOF) 
method [16], [17]. Although these algorithms are noniterative, 
their computation involves matrix inversions and a polyno- 
mial factoring or a solution of the generalized eigenvalue 
problem, which are still computationally inefficient compared 
to the proposed method. By taking into account the specific 
properties of the coefficient functions K± n which are real- 
valued and nonoscillatory, we can utilize more simple and 
efficient methods than the universal approximation techniques 
mentioned above; in particular, real exponentials are sufficient 
to approximate the coefficient functions and avoid any com- 
plex operations 1 . Thus, we have devised a simple relaxation 
algorithm based on curve fitting. The details of the procedure 
are given in Appendix A. Although this method is simple and 
iterative in nature, it converges to reasonable accuracy in a few 
iterations, and requires less computation time as compared to 
those for the previously mentioned methods. 

In general, the coefficient functions K “ n i have asymptotic 
values, and an analytical extraction of these values should 
be performed to increase the accuracy of the approximated 
functions or to reduce the computation time. The expressions 
for these asymptotic values can be easily obtained numerically 
or analytically. Furthermore, some of these coefficient func- 
tions are often zero, one, or symmetry of the others and these 
properties can also be explored to reduce the computation time. 

It can be seen that there exists a pole at 7 = 0 for a case with 
a presence of both top and bottom ground planes, and K± n i 
can no longer be approximated with exponential functions. 
To remove the singularity for this case, one can subtract 
the original Green’s function by the Green’s function for a 
homogenous case. Then, the remain part can be approximated 
as before, and for the homogeneous Green’s function a closed- 
form or a fast converging summation form can be used in the 
space domain [18]. The detailed discussion of a stripline case 
for the 2-D case is given in Appendix B. 

1 The complex-valued exponentials can be viewed as the complex images 

in the space domain in this paper. 


Green's Function 
In the Spectral Domain 



Fig. 4. Comparison of approximated and exact Green’s functions in the 
spectral domain. 

The Green’s function for the structure shown in Fig. 3 is 
approximated and compared with the exact Green’s function 
in the spectral domain. The dielectric constants and c 3 

are taken to be 5eo, and sq, accordingly. The thicknesses 
of the layers ti and *2 are 0.6 mm and 1.0 mm, and the 
observation and source points, y and y 0 » are 0.6 mm and 1.6 
mm, respectively. A maximum number of eleven exponentials 
were used to approximate the Green’s function. As shown in 
Fig. 4, the approximated results agree with the original one, 
and the maximum relative error was 0.0003. It is important to 
observe that although the exponential approximation might fail 
for the large argument case due to its fast decaying nature, by 
extracting the asymptotic value and the exponential factor from 
the coefficient function, the limiting behavior of the overall 
approximated Green’s function would still remain accurate. 

After approximating the Green’s function in the spectral 
domain, one can obtain the expressions for 2-D and 3-D spatial 
Green’s functions using the following identities: 


1 1 



The above identities can be easily derived by considering the 
potential due to the unit point or line charges 2 . Thus, in the 
space domain, the approximated Green’s function for 2-D and 
3-D cases, in general, can be written as 

G3D(r 1 ro) = jtt E /j 3D,± ( r 1 r °) < 15a > 

;= 1 

G 2D (HPo) = -^^/f' ± (p|p 0 ). 05b) 

2n£m 

Again, the superscripts + and — are used to denote the cases 
for y > Vo and y < y 0 , respectively. For j = 1 and y > yo , 

2 Equation (14a) can be viewed as the static version of Weyl’s identity. 
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2-D Green's Function 
In the Space Domain 



Fig. 5. Comparison of approximated and exact Green's functions in the space 
domain. 


the expressions for /^ D,+ (r | ro) and /J D + (r | ro) are given 
by (16a) and (16b), shown at the bottom of the page. Here, 
. denotes the asymptotic value of K+ n j and we have 
assumed that K+ nj is approximated by 


si 


^.»j(7)= £ Ciije---- 7 . > = 1-2,3, 4 (17) 


i=l 


where N+ nj is the number of exponential functions used 
to approximate K± nj . The approximated Green’s function is 
also compared in the space domain for the 2-D case. The same 
structure is used as in the previous one with £2 = £1 = 2eo> 
and y and y 0 were 0.6 mm and 1.6 mm, respectively. A 
maximum number of exponentials used in the approximation 
was five. The exact Green’s function is obtained by applying 
the image principle in the space domain. Fig. 5 shows the 
comparison. The maximum relative error was 0.0033. 

Finally, we note that considering the forms of (15) and 
(16) it can be seen that the exponentials used to approximate 
the Green’s function in the spectral domain correspond to the 
weighted images in the space domain. 


IV. Solution Method for the Integral Equation 

Once we obtain the approximated Green’s function, the 
integral (1) can now be solved by the moment method. First, 
discretizing the surface or the cross section of each conductor 
with polygonal-type elements, the unknown charge density can 


be expanded with the pulse basis functions over the discretized 
elements. Then, applying the point matching technique to each 
cell, (1) for 3-D problems becomes 


V, = 



G(xi,y„Zi | x,y.z)ds\ 


1 = 1,2 N t . 

(18) 


Here, qk is the unknown coefficient to be determined, and Nt 
is the total number of cells used to discretize conductors. Sk 
is the surface of the Jfcth cell, and point matching is applied at 
the center of the ith cell {x l ,y i ,z i ). The similar equation can 
be also obtained for 2-D problems. Assuming that all possible 
ground planes are held with the same potential and the ith cell 
is in the jth conductor, Vi is the voltage difference between the 
jth conductor and the ground plane; in the case of no ground 
plane, any one of the existing conductors can be chosen to be 
the ground reference to define the voltage. For the capacitance 
computation, V, takes the value of 1 or 0 depending on the 
excitation of the jth conductor. 

With manipulations in y, j/;, and the terms due to the 
exponential approximation, the surface integration in ( 1 8) can 
be put into the following forms: 



dV 

y/{x - x ') 2 + {y - y ') 2 + U - 2 ') 2 



Similarly, for 2-D problems we have 

J ln(i/( x _ x ') 2 + (l/ ~ Vf)dl' = J ln(|p - p\)dl' 

= f \n(P)dl'. 


(19b) 


The evaluations of (19a) and (19b) over an arbitrary polygonal 
patch and a line segment are well-known and the closed-form 
formulas are given in [19]. 

Now solving the system of linear equations obtained from 

(18) will give the unknown charge distribution in terms of 
its basis function. It is important to point out that without the 
closed- form expression of the Green's function, the integration 

(19) must be performed, theoretically, an infinite number of 


*3D,+ /_ i _ \ \ 

h ( 1 0j m ’ na y/(x - so ) 2 + (y + 2/0 - 2 d n y + (z ~z SJJ 

K, n.l . 

+ y* 1 

i = l ’ ’ y/(x- *o) 2 + (y + y o - 2 ~dn + i) 2 + (z - 2o) 2 

/i D,+ (p | Po) = 1 • In (\A X - *o) 2 + (y + Vo- 2d„) 2 ) 

jV + 

+ jr c^ n l • In (y / (i - X 0 ) 2 + (y + yo- 2 d„ + C,i) 2 ) • 


(16a) 


(16b) 
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TABLE I 


Computation 

Ddabare « al. [4] 


141.38 -21.493 -0.9924* 

-21.491 91951 -17.844 

-0.9023 -17.844 87.495 

(pF/m) 


' 14109 -21.765 -0.8920 

-21.733 93.529 -18.098 

-0.8900 -18.097 87.962 

(pF/m) 



Number of Basis Functions 

Fig. 6. Comparison of the present method with the spectral domain approach 
as a function of the number of basis functions. 

times. For most practical problems, the time to construct the 
matrix of the linear equation (moment matrix) takes the major 
portion of the computation time. We have significantly reduced 
this computation time by approximating the Green’s function 
with a finite number of exponentials. 

V. Numerical Examples 

A computer program was developed based on the above 
algorithm, and it is capable of handling an arbitrary number 
of dielectric layers and conductors and designed to read mesh 
data from a conventional mesh generator to allow computation 
of the complex geometries and meshes. The algorithm is tested 
for both 2-D and 3-D problems, and the equivalent crossover 
capacitance of two orthogonal microstrip lines is computed 
based on the integral equation with both 2-D and 3-D Green’s 
functions. 

A. 2-D Problems 

First, a simple microstrip capacitance is computed to verify 
the method, and the results are compared with those for 
the spectral domain approach (SDA) in Fig. 6. The spec- 
tral domain method used for comparison solves the integral 
equation iteratively in conjunction with the minimization in 
the boundary condition error. A more complex geometry, a 
three-conductor system in a layered medium, shown in Fig. 
7, is also computed. The number of basis functions used 
was 50 for each conductor. In both cases, the maximum 
number of exponentials used to approximate the coefficient 
function K± n t was 7. Comparison with the results in [4] is 
shown in Table I. In [4], the spectral-domain Green’s function 
is numerically integrated to obtain the corresponding space- 
domain Green's function using a Gaussian quadrature formula 
in conjunction with analytical asymptotic extraction. 



Fig. 7. Three conductors in a layered medium. All dimensions of conductors 
and spacings are identical. 



Fig. 8. Ten conductors in a layered medium. All dimensions of conductors 
and spacings are identical, and all units are in micrometer. 


TABLE II 


307.13 

-41.2 

-11.34 

-6.28 

-5.351 

•218.8 

■4.966 

-1.385 

-0.814 

-0.729 

-41.20 

319.6 

-28.04 

-7.79 

-4.987 

*5.025 

-216.95 

-3.54 

-0.985 

-0.666 

-11.34 

-28.04 

310.5 

-24.29 

-8.587 

-1.366 

-3.503 

-218.4 

-3.18 

-1.148 

-6.279 

-7.79 

-24.29 

303.5 

-24.73 

-0.798 

-0.946 

-3.164 

-219.4 

-3345 

-5.351 

-4.988 

-8.588 

-24.73 

290.5 

-0.708 

-0.627 

•1.12 

-3.325 

-221.3 

-218.8 

-5.01 

-1.363 

-0.796 

-0.709 

231.7 

-2.074 

-0.393 

-0.182 

•0.134 

-4.954 ■ 

-216.97 

-3.494 

-0.944 

-0.628 

-2.074 

232.0 

-1.19 

-0.255 

-0.135 

-1.382 

-3.532 

-218.4 

-3.157 

-1.12 

-0.393 

-1.19 

231.8 

-0.8752 

-0.242 

-0.813 

-0.982 

-3.174 

-219.5 

-3.32 

-0.182 

-0.255 

-0.875 

231.6 

-0.763 

-0.729 

-0.665 

-1.145 

-3.34 

-221.4 

-0.134 

-0.135 

-0.242 

-0.763 

231.3 


A ten-conductor transmission line system above a thick 
dielectric substrate, shown in Fig. 8, is also considered. The 
total number of 300 basis functions was used to represent 
the unknown charges, and nine exponentials were used to 
approximate each coefficient function of the Green’s function. 

Table II shows the computed results. The same structure is 
considered in [3] using the free-space Green’s function with 
the basis functions which incorporate the edge singularities 
of the charge near the comers of the conductors. In [3], data 
are obtained using 160 and 190 basis functions for conductors 
and dielectric interfaces, respectively. The methods used in 
[2] and [4] are also employed to compute the same structure 
in [3]. According to [3], the methods used in [2] and [4] 
resulted in nonphysical values, for instance, negative self- 
capacitance values. The method used in [4] took the CPU 
times of 89611.19 s on an IBM RS-6000 station with 300 
basis functions, whereas, the method in [3] took 458.67 s. On ) 

the same machine, our method took only 85.7 s of CPU time. , 

For the final 2-D example, a multilayered stripline case, J 

shown in Fig. 9, is also analyzed. Fifty basis functions are ! 

used to discretize each conductor. We have obtained the values S 
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Fig. 9. Two conductors in a layered medium with two ground planes. 
Dimensions of conductors are identical. 


k 




Fig. 10. Geometry of the test microstrip bend. 

Cn = C 22 = 217.07 pF/m and Cn = C 21 = -107.76 
pF/m. Data obtained from the Ansoft Maxwell software are 
C\\ — Cn = 217.65 pF/m and C 12 = C 21 - -108.24 pF/m. 


TABLE III 



Computation 

Silvester et al (20) 

Gupta et al [21] 

C„(fF) 

65.3 

72 

67.65 

C^rmalned 

<c„ nc^f • *)) 

0.359 

0.386 

0.363 


I I 
I I 


IK .A 

Ltzr 1 


, 1 i 1 

■ I j f | Conducior 2 
t 1 


Center Lines 


(a) 

Free Space 


frl 


(b) 

Fig. 11. (a) Planar and (b) cross-sectional views of the microstrip crossover. 


B. Equivalent Capacitance for a Microstrip Bend 

In this numerical example, we will use the 3-D Green's 
function to calculate the equivalent capacitances for a right- 
angle bend of a microstrip line. Although the equivalent 
capacitance or excess capacitance due to the right-angle bend 
discontinuities has been well-studied [20], [21], it is consid- 
ered here again to verify the present method. An accurate 
characterization of a microstrip bend involves the semi-infinite 
microstrip lines which, in turn, requires different expressions 
for the Green’s function. However, since the effects of the 
discontinuity are localized, it is expected that the equivalent 
capacitances can still be accurately calculated by truncating 
the semi-infinite lines with finite lengths. 

Referring to Fig. 10, the excess capacitance C ex can be 
defined by 

Cex = lim [Ct - C un [f(li + ^ 2 )] (20) 

'i.2— • 

where Ct is the total lumped capacitance and C un if is the 
capacitance per unit length for the uniform microstrip line. 
Setting V to 1 in (18), (20) can be directly written in terms 
of charges: 

C ex = Qex= Hm [Qt “ Qumfih + h)]- ( 21 ) 

Here, Qt is the total charge on the conductor, and Q un if is 
the charge per unit length for the uniform line. The excess 


capacitance was computed for w = 1.5 mm, h = 3 mm and e r = 
4.5, and l\^ and l[ 2 were 5 w and 20 w, respectively. A total 
of 456 nonuniform cells were used to mesh the geometry, and 
C U nif i which was calculated using the 2-D Green’s function, 
was found to be 62.12 pF/m. Table III compares the program 
results with those obtained from [20], [21]. The result from 
[21] contains a maximum error of 5%, while [20] contains 4% 
according to their analysis. The graph reading error must also 
be considered in [20]. 


C. Equivalent Capacitance for a Microstrip Crossover 

In the following, we seek the equivalent capacitances for 
a microstrip crossover. Fig. 11 shows the geometry of the 
crossover considered, and Fig. 12 shows the equivalent circuit 
representation. In Fig. 12, two lines are uncoupled except c m 
at the location of the crossover. Here, instead of constructing 
the integral equation in terms of the total charges, as we did 
for a microstrip bend case, we will formulate the integral 
equation in terms of the excess charges. The major drawback 
of the prior approach is that since the excess charge due to the 
discontinuities is usually much smaller than the total charge, 
the final accuracy in terms of the excess charge is much worse 
than the accuracy obtained for solving the total charge. 

Based on the equivalent circuit representation in Fig. 12, 
the integral equations after applying the boundary conditions 
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TABLE IV 


(22b) 


Here, the superscript T denotes the total quantity on the line 
while the superscript unif denotes the quantity per unit length 
of an isolated line. The integrations over the Green’s functions 
are symbolically written by {•,•), and the excitation voltage is 
assumed to be 1. The subscript ij used in the Green’s functions 
Gij denotes that the source is located at the ith conductor, 
and the boundary condition is applied at the jth conductor. 
The subscript ij used in the charge density functions q tJ 
denotes that qij is the charge density at the j>th conductor with 
excitation at the ith conductor. Now noting that (G t 2 t D , q£ nif ) 
is equal to V\ we can rewrite (22) as follows: 


(23a) 

(23b) 


Since the excess charges are expected to be localized, we 
can truncate the domain of the basis functions used to expand 
the excess charge density over distances /i and I 2 from the 
crossover junction (see Fig. 11(a)). Now using pulse basis 
functions and the point matching technique, (23) can be put 
into the following matrix form: 


r 0 1 


(^ 1 ? 1 + {^ 2 ? 1 ^ 12 ) 

,~( G i2 > 9 n if ) 

II 


~ (^ 2 ? > 922 lf ) J 

II 

( Gi ? W 2 l ) + ^ 2 ? ^ 2 ?) 

. + (@22 ' $22 m 


v n 

0 

0 

qjff ir 


M 1X 

m 12 ' 

q!i 

qS* 

0 

v„ 

0 unif 

L^n 

0 


[m 91 

m 22 

q'2 

qSL 



Erl = 2, and ^2 = 2 

£f j = 4. and £f 2 = 2 


69.17 (fF) 

59.748 (fF) 

4 

-55.40 (fF) 

-48.580 (fF) 

c{ 

-58.94 (fF) 

-53.968 (fF) 


Fig. 12. Equivalent circuit representation for the microstrip crossover shown 
in Fig. 11. 

on the surfaces of the conductors can be written as follows: 

(c®>qTi) + 0^2? ^ 12 ') 

.(^ 1 ? 'Q 11 / + ^ 2 ? » Qi: 2 / . 

. (g?P, iff") + (g;?, + (eg, ?;:) 

(^12* >921/ + {(*22 > 922 / . 

(^2?.922 if ) + (<???. 9l?) + ((?&?, flg)' 
(<?£> <7 2 u 2 nif ) + (G? 2 D , g? + (Gif, q%) _ 

( 


TABLE V 

Comparison of Results for Fig. 11 with i r \ — £ r2 =2 



Papatheodorou et al. 
[10] 

Computation 

c m overhic^ 

1.672 

1.776 

c\ 

-1.345 

-1.321 

c* 

-1.296* 

-1.513 

* This entry is incorrect, and it can be verified by considering 1 10, Fie. 4]. 


where 


rv»] y 

= - / g??(a 1 

p')dl' 


[ V 22]ij 

= - f G™(p, I 

J ci 

p')dV 

(24b) 

M k iij 

Qv 

CD 

II 

r')ds' 

(24c) 


<N 

8:? 

cr* 

8:? 

.JSLj 

II 

■ 3D,. 

/ 

q,“ nif 

— \n uni{ n unif 

- [y»,i 2 1 • • 

n Unif 

• * 9 i t , .Vj d . » 

]'. (24d) 


Here, [A] tJ? is the jth element of the ith row of the matrix A, 
and iV 3 D,i and A r 2 D,t are the number of patches and the number 
of segments used to discretize the surface and the cross section 
of the ith conductor, respectively. qjj nif can be obtained by 
solving the 2-D problems for each isolated microstrip with the 
excitation voltage set equal to 1. Then (24) can be solved for 
the excess charges. Finally the equivalent excess capacitance 
can be found by 


* V 3D,t 

c ii=Yl Area * • c* 

1 


(25) 


(24a) 


where Area* is the area of the fcth patch in the ith conductor. 

For a numerical example, the same structure used in [10] is 
considered, where h\ and /12 w cre 4 mm and 6 mm, the widths 
of both strips were 0.04 /ij, respectively, and I 2 were 10 
^ 2 » and 800 nonuniform patches and 40 pulse basis functions 
are used to solve the 2-D and 3-D problems, accordingly. 
Ttoo different dielectric configurations are considered and the 
results are shown in Table IV. The maximum number of 
exponentials used for each coefficient function was seven for 
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Fig. 13. Plots of the magnitudes of the excessive charge distributions for (a) q\\ , (b) (c) ^ (d) gif Crl = **• ^ Cr2 — x an< ^ 

y denote the transverse and the longitudinal, respectively. 


both cases. The result with e r \ = £ r 2 = 2 are compared 
with the results from [10] in Table V. In Table V, c\ and c' 2 
are the capacitances per unit length of the isolated wire lines 
with radii equal to 0.25 w and e r \ = £ r2 = 1* In [10], l\ 
and / 2 were 10 h\\ hence, their value of c m is expected to 
be smaller than ours as shown in Table V. The 3-D plots of 
excess charges are shown in Fig. 13. The Green’s function used 
in [10] is based on the image principle in the space domain 
and involves an infinite summation. According to [10], only 
ten terms were sufficient to evaluate the infinite summation 
for this particular structure, where lines are extremely narrow. 
However, for wider microstrip lines, the number of terms to 
compute the infinite summation will be much larger; this can 
be easily seen from the expression of the Green's function in 
[10]. 

VI. Conclusion and Future Work 

In the design of VLSI circuits, the capacitance computation 
plays an important role, and the present method can be 
conveniently combined with CAD tools. Since the Green’s 
function depends only on the processing parameters, such as 
the thickness and number of dielectric layers, the dielectric 
constants, and location of ground planes, the Green’s function 
approximation has to be performed only when the processing 
parameters are changed. Then, the approximation result can 
be used to compute capacitance parameters throughout the 
design period. 

In the computation of parasitic or excess capacitances, semi- 
infinite microstrip lines are often encountered, and to have 
accurate results, as mentioned earlier, the integral equation 
is often formulated in terms of the excessive charges which, 
in turn, requires the Green ’ s function due to the semi-infinite 
lines. The authors are currently studying the generalization of 
the present method for such a Green’s function in a layered 
medium. 

Appendix A 

The Real- Valued Exponential Approximation 
Based on the Relaxation of Curve Fitting 

First, we will assume that a function y(x) to be approxi- 
mated is real valued and nonoscillatory and, furthermore, that 
its asymptotic value is zero. The latter assumption can easily 
be satisfied if the function is limited at infinity. Our goal is to 


find the right-hand side of 

A’ N 

Viz) = f a {x) = ^2fi(x) = (A1) 

1=1 

Let w us for a moment assume that each first-order function f t 
approximates the original function y at some interval around 
one of the approximation points and is decreasing fast enough 
so that its value is negligibly small at the approximation points 
corresponding to larger values of the argument x. Then, we 
can safely determine one of the first-order functions, say fi, by 
neglecting contributions due to the other first-order functions, 
which are as yet unknowns to be determined. The parameters 
of f\ can be easily obtained by curve fitting two values of y 
for some large value of x. In a similar manner, we can find 
the parameters of the other first-order functions; however, this 
time we have to take into account contributions due to the 
previous ones which are already known. 

From the above argument, given 2N approximation points, 
the equations used to determine the parameters for the ith 
first-order function are then written by 

A _ 27rln(y,(2x2,-i)/y.(x 2 i)) (A2) 

£2;--£2j-1 

Ci = e- A ' x ’-‘ J/l (x 2 .-i) (A3) 

where 

yii x i) = yi-i( x j) - fi-ifa); j = (A4) 

In the above, yo(xj) is equal to y(xj). Let us now consider 

the case in which the value of a first-order function is not 
negligible at the other approximation points. In this case, if we 
perform the above procedure, there will be some difference 
between the original and the approximated function since 
we have ignored contributions due to some of the first-order 
functions. In such a case, to reduce this difference, we can 
iterate the above procedure including the contributions due to 
all other first-order functions which were obtained from the 
previous iteration. Thus, for the kth iteration, (A4) must be 
modified as 

Vi k) ( x j) = y( x i) - 53 fi k) ( x i ) ~ 13 

1=1 !=. + l 

j = 1 2JV. (A5) 
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v(j) » +e 1 
\1 + x 



x 


(a) 

Location of Pole vs. # of Iterations 



Number of Iterations 

(b) 

Fig. 14. (a) Comparison of approximated and exact values of (A6). (b) 

Convergence of the location of the smallest pole (-1) with the number of 
iterations. 


It can easily be shown that if y(x) had N distinct eigen- 
values, by iterating the above procedure, the approximated 
function will converge to the original function. However, in 
general, the approximating function will never be exact and 
the iteration must be stopped at some point at which the 
approximated function is optimal in the sense of curve fitting. 
Since our curve-fitting algorithm, most likely, has the biggest 
relative error for the values of x between £ 2 A r -i and xj.y* one 
can check the approximated value with the exact value at any 
point in this interval for the convergence criterion. In some 
cases, if the desired accuracy can not be achieved, then one 
should increase the number of exponentials. This case can be 
determined by checking the difference between the computed 
parameters of and 

The described method allows one to find the parameters of 
an approximating function one-by-one directly without solving 
a system of nonlinear or linear equations and does not require 
an original function to be monotonic. Moreover, the limiting 
values of the approximated function match exactly with the 
original function. 


Finally, to demonstrate the method, the following testing 
function is approximated and the results are shown in Fig. 14: 

y ( x ) = ~lf=r= + e_I - (A6) 

vl + x 

As illustrated in Fig. 14(b), to find the exact location of the 
pole at -1, we need a large number of iterations; however, 
since our goal is to approximate the overall function, only a 
few iterations were needed to approximate the function. To 
locate the poles exactly, one should use other methods such as 
those based on pencil of functions or the Prony approximation. 

Appendix B 


Strip Transmission Line Case 

Solving Poisson’s equation with separation of variables, 
the Green’s function for the strip transmission line case with 
the distance h between two ground planes can be written as 
follows (see [22]): 


C(x.» !✓.✓)= if Ian (12) 
7T£ ' n V h / 


n=l 

x sin 




V h 


(A7) 


The above expression for the Green’s function can be easily 
integrated analytically over the pulse basis function analyti- 
cally. This summation is quickly converging, and only a few 
first terms are needed; for instance, four digits of accuracy can 
be obtained using less than five terms for \x - x'\ jh > 0.5. 
However, when \x - x*\/h is small, the summation converges 
rather slowly, and the use of this Green’s function expression 
should be avoided. In such cases, we can use the following 
closed-form Green’s function: 


G(x.y | x\y') = - — In 
47T£ 


sinh 2 

Tr(i-x') 

+ sin 2 

^(y+y') 1 

- 

2 h 

2 h 


sinh 2 

Tr(x-i') 
2 h 

-f sin 2 

> 

1 ■* 



(A8) 


This expression for the Green’s function can be obtained using 
the conformal mapping and the method of images [18]. The 
major disadvantage of this expression is that the closed-form 
integration over the pulse basis function is unappealing; hence, 
the numerical integration scheme is unavoidable. Therefore, 
this form of the Green’s function must be only used when 
\x - x'\/h is small. As shown in [1], to integrate (A7) 
numerically, it is convenient to rewrite this expression by 
extracting the singularity as 


G{x.y | x'.y 1 ) = ^i-ln[(x - x') 2 + (y - y') 2 } 

+ g{x.y\ x',y') (A9) 


g(x,y | x',y') = —In 


[(x - x') 2 + {y- 

y') 2 ] 

Jsinh 2 

2h 

Ll' 

+ sin 2 

*(y+y')ll 

2 * JJ 

sinh 2 

2k 

4- sin 2 

w(y-y') 

2k 



(A 10) 
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where (A 10), shown at the bottom of the previous page. Now 
the first term in (A8) can be analytically integrated using 
(20(b)), and the second term can be integrated numerically. 
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Computation of the Equivalent Capacitance 
of a Via in a Multilayered Board Using 
the Closed-Form Green’s Function 

Kyung S. Oh. Jose E. Schutt-Aine, Raj Mittra. and Bu Wang 

Abstract — A method based on the quasi-static approximation for 
computing the equivalent capacitance of a via is presented in this paper. 
The geometry of a via consists of traces, pads and a perfectly conducting 
cylindrical rod; the via is buried in a multilayered dielectric medium with 
optional reference (ground) planes. The total number of traces, pads, and 
ground planes can be arbitrary, as well as the angles and cross sections. 
The method is based on the excess charge formulation of an integral 
equation applied in conjunction with the recently developed closed-form 
Green’s function. 

I. Introduction 

Although a via is one of the most common discontinuities en- 
countered in high-speed integrated circuits, u has not received as 
much attention as some of the other discontinuities, e g., open- 
end terminations, bends, and junctions. This is due mainly to the 
nonplanar and complex three-dimensional (3-D) geometry ot the via. 
which has often been simplified in the works published previously 
[|)_[3j. For instance, a via penetrating through a single reference 
(ground) plane with two wire traces has been considered in [11, 
and without any traces in [2], while a via above a reference plane 
with two wire traces but without a through-hole reference plane 
between these traces has been investigated in |3). A novel equivalent 
network model, which accounts for the frequency dependence, has 
been proposed in [4] and has also been applied to the problem 
of coupling between two adjacent vias in [5]. In [1] and [3|, an 
integral equation has been formulated in terms of the excess charge 
distribution to compute the equivalent (excess) capacitance. In this 
paper, this excess charge formulation is further generalized for vias 
with more complex geometries than has been analyzed hitherto and 
is applied in conjunction with the closed-form Green s function to 
analyze vias embedded in multilayered dielectric media. 

A closed-form Green's function for a multilayered dielectric 
medium was first introduced in [6). This Green's function utilizes 
a finite number of complex images and avoids the evaluation of a 
nested infinite senes expression required in the computation ot the 
exact Green's function for a layered medium. In [7], a closed-form 
Green s function based on weighted real images was proposed and 
was used to compute the equivalent circuit ot a stnp crossover. In 
[8], this closed-form Green's function based on weighted real images 
was further generalized to handle a semi-infinite line and then applied 
to compute various strip junction discontinuities. In this paper, the 
closed-form Green's function discussed in [8] is employed to compute 
the equivalent circuit for a via in a multilayered dielectric medium. 

II. General Statement of the Problem 

To illustrate the geometries of vias considered in this paper, a via 
comprised of three traces and one reference ground plane is shown in 
Fig. 1(a). In general, a via can pass through X g reference (ground) 
planes and X, traces, and X p pads can be attached to the via where 
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Fig. 1. (a) The geometry of a via. and (b) its surrounding multilayered 

medium. 

X g . Xt , and ,Y P are all arbitrary. The vias are embedded in a multi- 
layered medium consisting of X<i (arbitrary) dielectric layers, which 
can be backed by two optional reference planes as shown in Fig. 1(b). 
To distinguish between these optional reference planes and those 
associated with the via. we will reserve the term "reference plane" to 
designate an optional top or bottom reference plane, and use the term 
“reference conductor * to denote other reference planes. It is evident 
that the reference conductors must have perforations to avoid any 
contact with the vias; however, the two optional reference planes are 
assumed to be solid. To simplify the numerical computation, we as- 
sume that all of the conductors are infinitely thin and that the shapes of 
all the pads and perforations in the reference conductors are circular. 

A quasi-static equivalent circuit representation of the via shown in 
Fig. 1 is given in Fig. 2. This paper will only address the problem 
of computing the total equivalent capacitance CV . The method to 
compute the equivalent inductance of a via can be found in [9], 

In Section III, an integral equation is formulated in terms of the 
excess charge distribution using the closed-form Green’s function, 
and the method of moments (MoM) is subsequently employed to 
determine the unknown charge distribution. A detail discussion of 
the closed-form Green’s function and the corresponding expression 
can be found in [8]. In Section IV. several numerical examples are 
presented to verify the proposed method. 

III. Formulation of an Integral Equation 

An impressed potential on conductors results in free charge accu- 
mulation on the surfaces of conductors, and the electrostatic potential 
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TABLE 1 

Equivalent Capacitances for Vi as Shown in Fig 3. Units are in pF 



Fig. 3(a) 

Fig. 3(b) 


n 

£?=£ 0 


e f =4 5E n , e : =5Ae n 

Computed Result 

0.3841 

1.233 

9.952 

12.31 

Others 

0 3701 [3) 

1.28 [3] 

6 35(11 

7.85(11 


h 



Fig. 2. A equivalent circuit representation of the via shown in Fig. 1(a). 

o(r) at any point except inside conductors is then related to this 
surface charge density g(r) via the following integral equation 

o(r)= f G iD (r \ r' )q{r' )Hr‘ = (G 3D . q) (1) 

Ju 

where 12 denotes the surfaces of all conductors, including the via 
and reference conductors. G AU {r \ r ) is the 3-D closed-form Green’s 
function for a multilayered medium, and it accounts for polarization 
charges on dielectric interfaces and free charges on the surfaces of 
reference planes [8]. The integration is symbolically written as {-.*) 
to simplify the notation. Next, the charge distribution over the whole 
structure analyzed is decomposed into 7i(r).</p(r )•<//(/'). and q[ t { r ) 
where //,(r) is the charge density on the surface of a via hole, and 
q\ t { r ). q‘, ( /• ), and q u \r) are the charge densities on the /th pad, trace, 
and reference conductor, respectively. Equation (1) can be rewritten 
as 

A p 

o(r) = (G U) .<!,) + y^{G' l} .q' ft ) 

i 

+ ql) + qlf). (2) 

1=1 1=1 

Next, the charge density q\{r) is further decomposed into the uniform 
charge density '(r) and the excess charge density q\ xz ' ‘'(r) 

q't ( r ) = q)"" [ '( r) -F q, r ). (3) 

Here. 7| l,,lf,, (r) is the uniform charge density on the /th trace 
under the assumptions that it is infinite in both directions, no other 
traces are present, and the reference conductors have no perforations; 
this charge density is computed by solving an appropriate two- 
dimensional (2-D) problem. Since the reference conductors become 
uniform planes without any perforations for this 2-D problem, the 
potential distribution on the region above the reference conductor 
is not affected bv the region below it and vice versa. To solve 
for 7 ; nn ' 'in. it is then expedient to introduce a new medium 
surrounding the / th trace. As a consequence, the medium employed 


in the 2-D problem is generally different from that of the 3-D via 
problem, and could, in fact, be different for each trace. Once the 
appropriate medium has been chosen, q ) ni,i '(/*) can be obtained by 
using the method described in [7]. The resulting */J m,r,, (r) yields 
the capacitance per unit length for the /th transmission line in the 
equivalent circuit representation shown in Fig. 2. 

In the process of determining 7j ,nif '(/ ), the 2-D closed-form 
Green's function G 2U ‘{p \ po) is used to formulate an integral 
equation for a 2-D problem [7]. However, in the integral (2) for 
computing the equivalent capacitance problem, the uniform charge 
density q" ntL, ir) resides on the /th trace, which is only a semi- 
infinite line. It is therefore necessary to employ '<> | /<,.£) 

to compute the potential due to q" tili '(r) [8]. Using (3). (2) can be 
rewritten as 

•Vf 

, , \ "y/osemu unif.M 

O(r)- 2JG *<// ) 

< = i 

= (G 3 n .q,) + Y 2 (G 3D .^) 

r=l 

V 9 

+ <4>. (4) 

i=l r=l 

If we set the via potential to be Oo with respect to the reference 
conductors and planes, o(/ ) becomes Oo on the surfaces of the via 
hole. pads, and traces, and is equal to 0 on the surfaces of reference 
conductors. Hence, once q\"" [ ‘(r) has been determined, all of the 
quantities associated with the left-hand side of (4) can be considered 
to be known at the surface of the conductors, and the method 
of moments can be applied to solve (4). The various integrations 
appearing in (4) can be evaluated analytically for pulse-type basis 
functions using closed-form formulas given in (8). 

Once the unknown charge distributions have been determined, the 
equivalent (excess) capacitance C, can be obtained by using the 
following expression which involves the integrals of these charge 
distributions 

C. Oo = [ f V f q l v {r )dr 

Jiir fr * Jn p ., 

Af * f 

+ / q , ‘( r )dr* + / q' tf (r')th (5) 

where 12, is the surface of a via hole. S2 f , il ( l . and i 2, ; and are the 
surfaces of the /th pad. trace, and reference conductor, respectively 

IV. Numerical Examples 

We will now present two numerical examples to illustrate the 
application of the method presented above to the computation of the 
equivalent capacitances of two via structures, one with a reference 
plane and the other with a reference conductor. The detailed geome- 
tries of the two via structures are shown in Fig. 3. The computed 


IEEE 1 RAN SAC 1 IONS OS M!( ROWAVE THEORY AND TECHNIQUES VOL 44 NO 2. f KBKL ARY I'Wft 



(b) 

Fig. 3. (a) A two-irace via with a reterence plane and (b) a two-trace via 

with a reference conductor. All dimensions are in mm. 


excess capacitances for Fig. 3(a) with 1) fi — :j — -"0 an d 2) 
fl = 4 ; () and ii = for Fig. 3(b) with 3) r 1 = =2 = -Uo and 
4) fi = 4.5:i, and fj = 3.4?n are listed in Table I along with data 
obtained from (1) and [3]. In [1] and [3]. the strips were replaced 
by the equivalent wires of radii which are one-fourth of the widths 
of the strips. In our computation, the lengths of all traces have been 
truncated to 2.5 h, whereas the width of the reference conductor has 
been truncated to 1.5/, with h and I being the height of a via hole 
and the length of the traces. The truncation of traces and reference 
conductors is valid since the excess charge distribution decays rapidly 
as we move away from the center of a via. A total of 263 and 
687 unknowns were used for the vias shown in Figs. 3(a) and (b), 
respectively. As shown in Table I, the data for the via shown in 
Fig. 3(a) agree well with the published results. However, the data 
for the via shown in Fig. 3(b) are considerably different from the 
results reported elsewhere. Unfortunately, no experimental result for 
this structure is available to establish the relative accuracy of these 
results associated with Fig. 3(b). 

V. Conclusion 

A method to compute the equivalent capacitance of a via, which 
is based on an integral equation formulated in terms of the excess 
charge formulation, has been presented in this paper. The method is 


applicable to via geometries with or without ihrough-hole reference 
conductors. The recently developed closed-form Green s function was 
employed to circumvent the time-consuming evaluation of a nested 
infinite senes, required in the evaluation |3|. 
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Analysis of Edge Coupled Strip Inset Dielectric Guide 

Z. Fan and Y. M. M. Antar 


Abstract — The edge coupled strip inset dielectric guide is analyzed using 
the extended spectral domain approach. This structure, as compared to 
microstrip line, has several interesting features and can be very useful for 
microwave and millimeter wave applications. Validity of the approach is 
established by comparing numerical results with measured data. As many 
structural and material parameters can be chosen, a wide fundamental 
mode bandwidth and a broad range of characteristic impedances can be 
achieved, leading to great flexibility. The dispersion in fundamental mode 
propagation constants and impedances is found to be very low. With 
suitable choice of different permittivities for two dielectric layers, the 
same propagation constants for two fundamental modes can be obtained. 
This property is desirable for directional coupler applications. 

1. Introduction 

Microstrip line has been the most popular transmission medium 
used for constructing microwave and millimeter wave circuits [!]. It 
is well known that one of the problems with open microstrip circuits 
is the excitation of surface waves from discontinuities in the circuits 
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EFFICIENT TRANSIENT SIMULATION OF 
DISTRIBUTED INTERCONNECTS 
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I. Introduction 

The electrical performance of high-speed integrated circuits and digital networks 
strongly depends on the behavior of interconnects between various components of these 
systems. The prediction of such performance can only be achieved by the used of 
computer-aided design and simulation tools. The simulation of high-speed digital circuits 
has gained a significant role in the past few years since it is critical in the evaluation of noise 
levels, signal corruption and signal delay in fast switching circuits. This paper explores the 
various aspects and techniques for transmission line simulation; in particular, two different 
methods are described: the scattering parameter method and the optimal method. 

Two main aspects in the solution procedure must be distinguished, namely the line 
model and the transient simulation scheme. In the transient simulation step, the time- 
domain representation of the line model is obtained to yield the transient solution. 

Two different classes of models can be employed: circuit models and non-circuit 
models. Non circuit models cannot be directly implemented into a circuit simulator. Circuit 
models can be directly incorporated into a circuit simulator; they relate voltages and currents 
at the line terminals and are independent of the terminations. 

Further, two types of line characterization can be identified, namely closed-loop and 
open-loop characterizations. A closed-loop characterization such as Z- Y- H- and S- 
parameter characterizations results in complex oscillating transfer functions and transient 
characteristics. An open-loop characterization — in terms of the propagation functions — 
separates forward and backward waves and results in simple transfer functions and transient 
characteristics. 


IT. Scattering Parameters 

Scattering parameters have been extensively used in microwave applications for the 



measurement of high-frequency circuits as well as circuit analysis. The advantages of the 
scattering parameter approach reside in the flexibility in the choice of the reference system 
which can be adjusted to match certain characteristics of the network under study and 
simplify the analysis and the computations [l]-[3]. A direct consequence of this property is 
the greater computational stability and the simplicity of the closed-form solutions. More 
specifically, lossless and lossy transmission lines can be analyzed in conjunction with 
nonlinear terminations to yield solutions in which backward and forward waves are 
separated (see Fig. 1). 


Arbitrary Nonlinear 
Source Load 


Lossy 

and Dispersive 
Line 


Nonlinear Arbitrary 
Load Source 



Port 1 


□— z 2 (t) — On 

x=l V 

Port 2 


x=0 


Fig. 1. 


Consider a transmission line of length 1 with complex characteristic impedance Z c and 
propagation velocity v. This transmission line is connected to two ideal lossless reference 
lines as shown in Fig. 2, it can then be described in the frequency domain through its 
scattering parameters using 


Bi = Si i A i +S12A2 (l a ) 

B 2 = S 21 A 1 +S 22 A 2 (lb) 


where Aj and A 2 are the incident voltage waves in the reference line of characteristic 
impedance Z 0 . B 1 and B 2 are the reflected waves due to A 1 and A 2 respectively. It can be 
shown that 


(1 . e -2ja>i/v)p 

11 22 1 . e -2jcol/vp2 


S 12 = S21 = 


(1 - p 2 )e-i^ v 

1 . e -2jcol/vp 2 


(2) 


where 


P = 


Zc ~ Zp 

Zc + Zq 


(3) 




In the time domain this correlates to (lowercase characters for time-domain variables) 


b] (t) = s | ,(t)*a i(t) + s i 2 (t)*a 2 (t) (4a) 

b 2 (t) = s 2 i(t)*a i(t) + S22(t)*a 2 (t) (4b) 

where Sjj(t) is the time-domain Green's function associated with Sjj. These equations are 
combined with the terminal conditions 


ai(t) = T i(t)gj(t) + n(t)b,(t) 


(5a) 


a 2 (t) = T 2 (t)g 2 (t) + r 2 (t)b 2 (t) 


(5b) 
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Fig. 2. 


where Tj(t), T 2 (t) are the transmission coefficients, T i(t) and T 2 (t) are the reflection 
coefficients and gj(t) and g 2 (t) are the time-domain expressions for the voltage sources. 
Tj(t), T 2 (t), Ti(t) and T 2 (t) are time-dependent and are related to the terminations Z\(l) and 
Z 2 (t) of Fig. 1 which may be the time-domain representations of reactive elements or 
linearized active circuits. The convolution operation is defined as 

t 

Sjj(t)*aj(t) = Js ij (t-T)*a j (T)dT 

o (6) 

When time is discretized the convolution becomes 

Sij (t)*aj(t)= is ij (t-t)a J (T)Ax (7) 

T=1 

where Ax is the time step, 
we isolate the term containing aj(t) 




(8) 


t-i 


Sjj(t) * aj(t) = s ij (0)a j (t) + I Sjj (t - x)aj(x)Ax 

T=1 


we then define the history term as 


Hjj(t) = Xsij(t-T)aj(T)Ax : History 

T=1 


Define s'jj(O) =sjj(0)Ax we finally obtain 


(9) 


bi(t) = s' i i(0)aj(t) + s'i2(0)a2(t) + Hn(t) + H]2(t) (10a) 

b2(t) = s'2i(0)ai(t) + s’ 2 2(0)a 2 (t) + H2i(t) + H22W (10b) 

When combined with the termination conditions of (5), the solution for the independent 
waves can be extracted to yield 

ri-r2(t)s' 2 2(0)][Ti(t)gi(t)+ri(t)Mi(t)] , r,(t)s'i2(0)[T2(t)g2(t)+r 2 (t)M2(t)L Mla> 

ai(t)= A(t) + A(t) (Ua) 


a 2 (t) = 


n-ri(t)s'n(0)] T2(t)g2(t)+r 2 (t)M 2 (t)] ^(Os^ia^rriWgKo+rKQMKt)] 

A(t) A(t) 


A(t ) = [i-r,(t)s' n (0)] [i-r 2 (t)s'22(0)]-ri(t)s'i2(0)r 2 (t)s , 2i(0) (lie) 

with Ml (t) = Hi 1 (t) + H12 (t) and M 2 (t) = H 2 i (t) + H 22 (t). The total voltages and 
currents at the ports can be recovered from the voltage waves using the relations 


vi(t) = ai(t) + bi(t) 

(12a) 

v 2 (t) = a 2 (t) + b 2 (t) 

(12b) 

a,(t) b,(t) 

1|<t) - z„ • z„ 

(12c) 

1 

II 

4—1 

(12d) 


In general, the frequency domain scattering parameters can be calculated from the line 



parameters using equations (2); next an inverse Fourier transform permits to extract the 
time-domain Green' s functions which can then be used to calculate the voltage waves. 
Most of the computations reside in the evaluation of the convolution history terms as 
expressed by (9). 


TIT. The Optimal Method 

In [4] the problem of distributed line simulation was analyzed to develop the optimal 
method, resulting in the maximum efficiency, accuracy and practical applicability with 
respect to the transient simulation of digital circuits. The method handles, in the same 
straightforward manner, uniform and nonuniform multiconductor lines with constant and 
frequency-depend parameters, including real lines characterized with frequency- or time- 
domain data samples. Moreover, an excellent simulation accuracy in the full frequency/time 
range, from zero to infinity, is attained with the minimal number of the data samples. The 
resulting line model can be directly used in a circuit simulator; the method supports variable 
time stepping and has linear computational complexity. The efficiency of the optimal 
method allows for an accurate simulation of realistic circuits, containing thousands of 
multiconductor nonuniform frequency-depend lines and nonlinear active devices, with 
virtually no increase in the simulation time compared to the simple replacement of 
interconnects with lumped resistors. 

As shown in [4] and [5], to achieve the maximum efficiency, accuracy and practical 
applicability, transient simulation of distributed lines should be based on: 

• a time-only formulation; 

• the open-loop characterization; 

• a device model which does not require to introduce current variables; 

• indirect numerical integration; 

• the difference approximation based on interpolation. 

The detailed description of the optimal method and its application to uniform and 
nonuniform coupled lines can be found in [5]-[9]. Here, we shall confine ourselves to a 
single uniform line case. 

For single uniform lines, the open-loop element characteristic is given by 


Il=Y c ((0)Vi-Gi, 
I2 = Yc(co)V 2 -G2, 


(13a) 

(13b) 



ii h 



Fig. 4. A circuit-diagram interpretation of the line element characteristic. 


where 

Gi=e'Y*[2l2 + G2], (14a) 

G 2 = e-Y 1 [2Ii+Gi], (14b) 

Y c =1/Z C is the complex characteristic admittance of the line, and e*^ is the propagation 
function. Conventions for the terminal voltages and currents are shown in Fig. 3. Figure 4 
shows a circuit-diagram interpretation of the element characteristic (13). 

The difference approximation [6] is a general method for applying numerical integration 
to systems whose differential equations are not available (such as systems characterized 
with discrete samples of their time- or frequency-domain responses, or transcendental 
transfer functions). The first step of the difference approximation consists of approximating 
the responses with the corresponding responses of a system with known differential 
equations. For instance, for the scalar frequency-domain difference approximation in the 
parallel canonic form, we express the propagation function in terms of rational function as 





where x is the propagation delay, A] is the final value, the aj’s and G) c j’s are the cut-off 
frequencies and partial-expansion coefficients respectively. We also approximate Yc(«) as 


Y c (co) = 



a 2i 

l + j(D/CO c2i _' 


(16) 


Note that there is no delay term since the transfer relationship of the impedance is 
instantaneous. 

We use interpolation-based methods [6], [9] to perform the approximation. One of the 
methods — the relaxation interpolation method in the parallel form — approximates the 
original function in terms of a sum of basis functions 

L 

g(©)-5>i«o). (17) 

i=l 

where g(co) is the function to be approximated and <J)j (oo) are the basis functions. It is 
assumed that the finite values of the basis functions and the finite value the original function 
are zero. Let each basis function approximate the original function at some interval around 
one of the approximation points and be decreasing so fast that its value is negligibly small at 
the approximation points corresponding to larger values of the argument. For the first basis 
function this assumes 


^(co^sgfco,). 


(18) 


Now, we can write for the next basis function 


4> 2 (co 2 ) = g(co 2 ) — <J)i(co 2 ), 

and more generally 

i-l 

<J>i (CDj ) = g(COj ) — 4>j(©i). 

j=l 


(19) 


( 20 ) 


After the first pass, we can repeat the above procedure, taking into account all the basis 
functions obtained from the first pass: 



( 21 ) 


L 

(JjjCcOj ) = gCcOj ) - 21 <t>j(cOi). 


j=l 

j*' 


After this second pass, the difference between the original and approximating functions 
will be reduced, because an influence of all the terms was accounted more accurately during 
the calculation of the approximation parameters. The procedure may be repeated to achieve 
a desired accuracy. 

The following summarizes the steps involved in applying the real-part parallel relaxation 
interpolation method to Y c (to). In this case, Y c ((0) is being approximated as 

a ; 


Let g be the real part of Y c , then 


L 

Y c (to) = | A +2 ■ 

it? l + J<o/®d 


L 

g(CO) = | A + £ 

i=l 


l + (o 2 /a)cj 


1. Choose approximation interval: 2L+ 1 frequency points. 


2. Sample original function at the frequencies 0) t = G) n 


2 L 


3. Reverse the order of the samples. 

4. Initialize ai = 0, coi = 1 , i = 1,2, . . . , L. 
REPEAT until a desired accuracy tolerance is attained 


( 22 ) 


(23) 


, k = 0, 1, ..., 2L. 


A = g(CD r 


L 

.)-! 




j=l 1 ® max / ®cj 


(24) 


FOR 

i runs as 1,2, ...., L 

w 2 _ Qj (° 3 4 2i ) CO li ~ Qj( C0 2i-l) C0 2i-1 

0i((O 2 i-i)-0i(®2i) , (25) 


a; = 0i(CO 2j ) 


1 + 


«2i 

°>d J 


(26) 


where 
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Fig. 5. Comparison of the original (thin continuous curves) and approximate (thick dashed 

curves) functions. 


L 

0i((O k ) = g(0 k ) — A-£ 


j=i 

j*i 


1 +C0 k /(Ocj 


(27) 


End of FOR loop, 

End of REPEAT loop. 


Figure 5 shows a typical 5th-order approximation of the dynamic part of the propagation 
function. 

The indirect numerical integration [7] is a class of numerical integration methods based 
on the analytical relationship between the continuous and discrete domains. It has ideal 
accuracy, convergence and stability properties. The approximate rational polynomial 
transfer functions represent systems described with finite-order ordinary differential 
equations, and can be discretized, and simulated recursively in the time-domain using 
numerical integration. For instance, for a system with the transfer function of the form 


A + 


y ^ 

h, l+iw 


£i 

j(D/C0 ci 


Jan 


H(f) = 


(28) 





the step-response invariant indirect numerical integration formulas in the parallel canonic 
form are given by 

L 

y[nT] = Ax[(n - k)T] + Xy pi [nT], (29) 

i=l 


where 


y pi [nT] = a s (l- e-^ T )x[(n - k - 1)T] + e^^y^ [(n - 1)T], (30) 

x is the excitation, y is the system response, T is the time step and k is the discrete-time 
delay. 

After the characteristic admittances and propagation functions have been discretized, the 
frequency-domain element characteristic (13) turns into the discrete-time element 
characteristic — the companion model — which is used by a circuit simulator to represent 
distributed lines, and is recursively updated at each time iteration via numerical integration 
formulas. 

Figure 6 presents the comparison of the simulation results obtained with the optimal 
method and segmentation method (separating delays and losses). The network and line 
parameters are: Rl=50 Q, R2=l k£2, C=39 pF/m, L=539 nH/m, R=125 QJm, G=0, 1=0.675 
m. The following voltage pulse with linear rise and fall was taken for the excitation: El=4 
V (start time t s =5 ns, rise time t r =l ns, pulse width t w =20 ns, fall time t,=l ns), E2=0. The 
time-step T=0.5 ns. As can be observed, the results are in excellent agreement. 

IV. Conclusions 

Two methods of distributed line simulation were discussed — the scattering-parameter 
method and the optimal method. The scattering-parameter method permits a straightforward 
implementation of voltage wave solutions but is computationally expensive as a result of 
convolution calculations. The optimal method results in the maximum efficiency, accuracy 
and practical applicability with respect to the transient analysis of digital circuits. The 
optimal line model can be directly used in a circuit simulator, and lines can be characterized 
with frequency or time-domain data samples. The efficiency of the optimal method allows 
for an accurate simulation of real circuits, containing thousands of multiconductor 
nonuniform frequency-depend lines and nonlinear active devices, with virtually no increase 
in the simulation time compared to the replacement of interconnects with lumped resistors. 




E2 
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Fig. 6. Network configuration and comparison of the simulated waveforms obtained using the 
segmentation method (thin continuous curves) and the optimal method (thick dashed curves). 
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Abstract — An optimal technique (In terms of computational 
efficiency, accuracy, and practical applicability) for the tran- 
sient simulation of distributed multiconductor RLGC lines was 
developed by applying the indirect numerical integration 
method (the most efficient and accurate transient simulation 
method) to the device model for a distributed line (the most ef- 
ficient and practically applicable model). The technique can be 
directly used In a circuit simulator. The computational 
complexity is linearly proportional to the number of time-steps. 
Nonuniform lines, lines with frequency-dependent parameters, 
and real lines, characterized with a few samples of time- or fre- 
quency-domain measurements or electromagnetic simulations, 
can be accurately and efficiently simulated. 

I. Introduction 

The transient simulation of distributed lines has gained 
special importance with the development of high-speed digi- 
tal electronics. Previously developed techniques are quite di- 
verse, but not always computationally efficient and applicable 
to real circuits in which distributed lines are surrounded by 
thousands of active devices. In [1] an attempt was made to 
develop a technique that would be optimal in computational 
efficiency, accuracy, and applicability to real problems. The 
results obtained in [1] are summarized in this paper. 

II. System Model for a Distributed Line 

The system-diagram representation of a single line is 
shown in Fig. 1. W Vf and W vb represent the forward and 
backward voltage propagation functions, respectively; Tj, T 2 
and r lf T 2 are the near- and far-end transmission and reflec- 
tion coefficients, respectively. 

The system-diagram representation of a multiconductor 
line is given in Fig. 2. Wvmf an ^ ^Vmb m diagonal 
forward and backward modal voltage propagation function 
matrices; T\ t T 2 and Tj, F 2 represent the near- and far-end 
transmission and reflection coefficient matrices; and Mvf» 
Mvb are the forward and backward modal voltage transfor- 



Fig. I . System model for i single line. 
0-7803- 1254-6/93S03.00 © 1993 IEEE 
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mation matrices. 

The models reproduce the general relationship between the 
physical processes of wave propagation, reflection, and 
coupling in a distributed system and can represent arbitrary 
distributed systems such as uniform and nonuniform trans- 
mission lines, waveguides, and plane-wave propagation. For 
distributed RLGC lines, the mathematical justification of the 
models was obtained [1]. The models can be used for system 
analysis of distributed networks. They also allow one to ob- 
tain the solution for a distributed system without complex 
mathematical derivations. 

From Figs. 1 and 2, it can be observed that a distributed 
line along with the terminations forms a feedback system. 
This explains why a global-parameter characterization (such 
as Z-, Y-, //-, or 5-parameter characterization) results in com- 
plex, non-monotonous functions for the parameters. A direct 
characterization in terms of open-loop transfer functions (the 
propagation functions) leads to simpler transfer functions and 
transient characteristics and should be used for the transient 
simulation. 

m. Device Model for a Distributed Line 

Distributed line models can be divided into circuit and 
non-circuit models. Non-circuit models cannot be directly 
inserted into a circuit simulator, and, consequently, can not be 
applied with efficiency for the transient analysis of real cir- 
cuits. The system models described in the preceding section 
may serve as an example of non-circuit models. 

Circuit models can be directly placed into a circuit simula- 
tor, and, therefore, are of prime practical interest. They relate 
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tion. The resulting difference equation is called the differ- 
ence model of the continuous system [4], Depending on the 
discretization technique, one recognizes the direct and indi- 
rect numerical integration methods. 

For the direct methods, the differential equation is dis- 
cretized using approximations for integrals or derivatives. 
Indirect numerical integration methods are based on the ana- 
lytical relationship between discrete samples of the system 
excitation and response (analytical relationship between the 
continuous and discrete time domains). The relationship is 
obtained using the time-response invariance method. For this 
method the excitation is assumed to be piecewise of a certain 
form, for instance, staircase, piecewise linear, piecewise 
quadratic, cubic, etc. Then, the general solution of the system 
differential equation can be found and expressed in terms of 
the discrete samples of the excitation and response (by using 
the samples as boundary conditions). The resulting differ- 
ence model is then used to recursively compute the response 
in the time domain [4]. 

The indirect numerical integration provides the exact sys- 
tem response to an excitation which is piecewise of the same 
form as the difference model invariance. For instance, the 
ramp-invariant difference model provides the exact system 
response to a piecewise linear excitation such as a pulse with 
linear rise and fall. The indirect numerical integration has 
ideal stability and convergence properties, because, in the ab- 
sence of an excitation, it always provides the exact system re- 
sponse, irrespective of the value of time step and form of in- 
variance. It is also more accurate than the direct method. 

From the foregoing discussion one can conclude that an 
optimal technique (in terms of computational efficiency, ac- 
curacy, and practical applicability) for transient analysis of 
distributed lines should be based on the device model and the 
indirect numerical integration. This combination was first 
proposed by A. J. Gruodis and C. S. Chang [5], but using the 
direct numerical integration method. 

Due to the distributed nature, modules of a transmission 
line model can not be exactly described with finite-order or- 
dinary differential equations. Consequently, an approximate 
ordinary differential equation representation has to be found 
before numerical integration can be applied. This can be 
achieved by approximating a frequency- or time-domain re- 
sponse with the corresponding response of a system described 
with ordinary differential equations. The following differ- 
ence approximation procedure was developed that allows one 
to apply indirect numerical integration to an arbitrary system. 
In general, a system may contain an attenuation component 
A, a delay component t, and a dynamic part. The difference 
approximation procedure consists of one step: represent the 
transfer function in the form 


W(jU3) = 


A + f — 

~ l + ;Ci)/0) c( 




or the transient characteristic in the form 


(3) 


h(t) = 






u(r-t) . 


J 


(4) 


Then, the step-invariant difference model is given by 

L 

>•(/„) = Ax(t n -t) + £>•,('«) . (5a) 

i-l 

where 

y,(r„) = fl 1 (l-e-° , « T »)x(r„_, -t) + e-“« r " (5b) 

and the ramp-invariant difference model is 

L 

y(r„ ) = Bx(t n - x) - (r„ ) . (6a) 

1 = 1 

where 

^•(l-e - ®"' 7 ") 

yiOn ) = „ r 4*('n " t) “ *(/„-! - X)] 

+ e~ Wc ' Ta y, (f n _i) • (6b) 

A and B are the final and initial values of the transfer func- 
tion, x(t) and yit) are the excitation and response, TV^iHn-i is 
the value of the time-step at n-th iteration, and u(r) is the unit- 
step function. Further increase in the order of the time-re- 
sponse invariance results in a difference model of higher or- 
der than the system differential equation without a significant 
increase in the difference model accuracy. Eqs. (3H6) cor- 
respond to the parallel canonic form of the difference model. 
Other forms (including direct and cascade canonic ones) can 
be found in [1]. 

In general, second-order terms representing complex-con- 
jugate poles of the transfer function have to be added to the 
sum in (3) and the corresponding harmonic terms to the sum 
in (4). These terms were omitted here, because the elements 
of the device model for a distributed line represent aperiodic 
systems, and, therefore, have to be approximated using only 
real poles. 

The frequency-domain difference approximation procedure 
is more general, because it can directly handle lines with arbi- 
trary frequency-dependent parameters or lines characterized 
with frequency-domain measured data. The time-domain dif- 
ference approximation procedure should be employed only if 
transient characteristics are available. For a single RLGC 
line, the analytical expressions were obtained for the transient 
characteristics and limiting values (A and f?) for all the 
modules of the system and device models [1]. 

The difference approximation procedure is applied to the 
characteristic admittances and propagation functions. In the 
case of a multiconductor line, all entries of the characteristic 
admittance and propagation function matrices need to be ap- 
proximated. The resulting time-domain device models have 
the same form as the frequency-domain models shown in 
Figs. 3 and 4, At each time iteration the instantaneous values 
of the conductances and current sources are recursively com- 
puted from the old values using the corresponding difference 
models. If the step-invariant difference model (5) is em- 
ployed, the conductances are constant throughout the simula- 
tion (because the first-order difference models (5b) do not 
contain the current value of the excitation), and only the cur- 
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the SIBC, it is represented by a convolution integral in the 
time domain. Tesche [6] has formulated a time-domain integral 
equation based on this convolution integral and pointed out 
that the direct computation of this convolution integral is 
impractical due to the large computation time and storage 
requirements. On the other hand, Maloney and Smith [7] have 
applied the SIBC to the finite-difference time-domain method. 

To overcome the computational difficulties associated with the 
convolution integral, they use an approximate recursive for- 
mula. Their implementation, however, requires the exponential 
approximation to be performed prior to the FDTD simulation. 
Beggs et al. [81 have eliminated this preprocessing uine by 
performing the exponential approximation based on the high 
conductivity surface impedance approximation. It should be 
noted that when the SIBC is applied to the finite-difference 
time-domain (FDTD) method, in addition to the reduenon of 
discretization space, further computational saving is achieved 
due to the larger discretization step compared to that for the 
nonreduced original problem. A detailed discussion of this 
computational saving is given in [7] and [8J 

In this paper, the implementation of the SIBC in the FDTD 
method is considered for a lossy dielectric half-space and 
a thin lossy dielectric medium. For a lossy dielectric h 
space, the normalized impedance function for a lossy medium 
is approximated in the frequency domain using a senes of 
first-order rational functions. Since the normalized impedance 
function used in this paper is independent of medium proper- 
ties, the rational approximation has to be performed only once 
The results of this approximation are tabulated for rauonal 
functions of various orders. Then, using the approximate 
normalized impedance function and assuming the waves to be 
piecewise linear in time, a closed-form expression is derived 
to evaluate the time-domain convolution integral recursively. 
Thus, the presented formulation is numerically more efficient 
than that in [7] because the preprocessing time for the expo- 
nential approximation is eliminated and more accurate than 
that in [8) since the high conductivity approximation is not 

“^Although several methods have been proposed to model a 
thin lossless dielectric shell for the FDTD method PH 12], 
little work has been performed for a perfect electric conductor- 
(PEC-)backed P lossy dielectric shell in spite of its frequent 
occurrence in practical problems. A full time-domain imple- 
mentation of the SIBC, which accounts for the frequency 
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Computation of Excess Capacitances of Various Strip 
Discontinuities Using Closed-Form Green’s Functions 

Kyung S. Oh, Jose E. Schutt-Aine, and Raj Mittra 

Abstract — An efficient quasi-static method to compute excess (equivalent) 
capacitances of various strip discontinuities in a multilayered dielectric medium is 
presented. The excess charge distribution on the surface of a conductor is 
obtained by solving an integral equation in conjunction with closed-form Green’s 
functions. A complete list of expressions of the closed-form Green’s functions 
for a point charge, a line charge, and a semi-infinite line charge is presented. An 
open end, a bend, a step junction and a T junction are considered as numerical 
examples. 


I. INTRODUCTION 

Q uasi-static analysis is often performed to characterize strip discontinuities when the 
dimensions of the discontinuities are much smaller than the wavelength. Under the quasi- 
static analysis, the dominant effect of strip discontinuities is fringing fields due to the physical 
irregularities of discontinuity geometries. The modeling of these fringing fields in terms of an 
excess capacitance is discussed in this paper. 

Numerous papers have been published to compute excess capacitances of various microstrip 
discontinuities, and a summary of popular methods can be found in [1]. The most successful 
approach is one based on the formulation of an integral equation in terms of the excess charge 
distribution, which was first proposed by Silvester and Benedek [2] and has been applied to 
analyze various microstrip discontinuities [2]-[5]. The Green’s function for a layered medium is 
employed in this approach. For N dielectric layers, the expression for this Green’s function would 
consist of an N-\ nested infinite series [6]; hence, in practice, this form of the Green’s function 
may not be applied to a multilayered medium. Recently, Sarkar et al. [7] solved discontinuity 
problems for a multilayered medium using the free-space Green’s function, but additional 
unknown charges (over unknown charges on the surface of a conductor) had to be placed on the 
dielectric interfaces and the top ground plane to model the polarization charge and the free charge. 
Although such inclusion of unknowns may be tolerable for 2-D problems, it is computationally too 
burdensome for 3-D problems. 

In this paper, the closed-form Green’s function discussed in [8] is employed to formulate an 
integral equation in terms of the excess charge distribution. A complete list of expressions of the 
closed-form Green’s functions for a point charge, a line charge, and a semi-infinite line charge 
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with or without a top ground plane is presented in this paper. The presented method requires 
neither additional unknowns to model dielectric interfaces and the top ground plane nor evaluations 
of any infinite series except for cases where the top ground plane is present. When the top ground 
plane is present, using the closed-form Green’s function is still numerically advantageous since the 
nested infinite series in the expression of the usual Green’s function becomes a simple infinite 
series without nesting. 

II. CLOSED-FORM GREEN’S FUNCTIONS 


Closed-form expressions of the electrostatic Green’s functions for a point charge, a line 
charge, and a semi-infinite line charge are derived in this section based on the approach used in [8], 
Consider N dielectric layers which are backed by a ground plane as depicted in Fig. 1. The Mh 
layer is either a half-space or terminated by an optional top ground plane. All dielectric layers and 
ground planes are assumed to be planar and infinite in the xz-plane. The electrostatic Green’s 
function in the spectral domain is described by the following closed-form formula [8]: 


G(y,y\r 0 ) = —^—(Kf(y,m,n)e y(y+y <’ 2dn) + K\ (y,m,n)e y{y y«+W m _ , d n )) 


+Kl{y,m,n)er ( - y+yo) + Kt{Y,m,n)e^-y-yo +2d ^) 


y^y 0 ( la > 


G(y,y\r 0 ) 


l 

2 £ m r 


+K 3 (y,m,n 


■( Kf(y,m,n)e Yiy+y ° 2dm) + K 2 (y,m,n)e Y{y y<>) 
yy(-y+y 0 + ^ d n-\~ d m)) + K^ (y,m,n)e Y( ~ y ~ y ° +2dn - l)] ) 


y^y a ( lb > 


where G is the spectral-domain Green’s function, and r and r 0 are the observation and source 

points located in the nth and mth layers, respectively. The superscripts + and - are used to denote 
the cases for y>y 0 and y<y 0 . The expressions for the four coefficient functions Kf can be 

found in [8]. A closed-form expression of the Green’s function in the space domain is obtained by 
approximating these four coefficient functions Kf using exponential functions. It is important to 
mention that although Kf is dependent on m and «, it is not a function of y and y 0 ; hence, the 
approximation can be performed without any prior knowledge of the junction geometry. The 
following two subsections detail the derivation of the closed-form Green’s functions due to a point 
charge, a line charge, and a semi-infinite line charge with or without the top ground plane. 

A. Closed- Form Green ’s Functions for Geometries Without the Top Ground Plane 
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When there is no top ground plane, the four coefficient functions K.f in (la) or (lb) are 
nonoscillatory and smooth functions of y ; hence, each coefficient function can be sufficiently 
approximated with real- valued exponential functions as follows [8]: 


Kf{m,n,y) 


N„ 


C ±J e 


+ / 

a * Y 

■ * 


7=1 


i = 1,2, 3, 4 (2) 


where N ^ n , denotes the number of exponential functions used in the approximation of K, , 
which typically ranges from 5 to 10. The exponential functions used in the spectral-domain 
approximation can be physically interpreted as weighted images in the space domain [8], 
Compared to the exact Green’s function in the space domain, which consists of an infinite number 
of images, the expression for the closed-form Green’s function consists of only a finite number of 
weighted images. 


Once the exponential approximation is performed in the spectral domain, the closed-form 
Green’s functions for 2-D and 3-D can be obtained in the space domain by using the inverse 
Fourier transformation formulas (14a) and (14b) in [8], and the resulting expressions are given by 

0 1D (p\p 0 )=-^-bf D ^P\Po) W 

/ = j 

<»> 

/ = i 

For t'=l and y > y 0 , the expressions of j ? D,± and f^ D,± are given by 


i 


fl D ' + {p\po)= X C m,i,l • ln N<* - x o) 2 +(y + y 0 - 2d n + “mil f 
7=1 


N, 


r+J 

u m,n,l 


f 3D,+ ( J r)= y 

1 0 j=l y/(x - X 0 ) 2 + (y + y 0 - 2 d n + a^) 2 ZoY 

Similar expressions can be obtained for ^ 2D,± and for other values of i. 


(4a) 


(4b) 


To derive the Green’s function for a semi-infinite uniform line charge, the auxiliary Green’s 
function for a line charge with polarity reversal is employed [2], Consider a uniform line charge. 
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which starts from z = % and is infinitely extended in the positive z-direction; then the Green s 
function for a semi-infinite line charge G sem can be expressed as 

G Itmi (r |r„,£)= i[c 2D (p|p„)+ G”(r\r c 4)\ (5) 


where G semi is the Green’s function for a line charge with an abrupt polarity reversed from minus 
to plus at z = |. The expression for G p is obtained by integrating the potential due to a point 

charge [2]: 


C p (^) = -jG 3D (r]r,,)+jG 
-°° £ 


3D 




( 6 ) 


m /=1 


Again, for i= 1 and y > y 0 , f p ' is given by 


<.„,i 

fr + (r\r o ,0= I C^ ln 

j = 1 


(x-x 0 ) 2 +{y + y 0 -2d n +a^ J JlA f +(z-Q +(z~g) 
(x - * 0 ) 2 + (y + y„ - 2d„ + fl^ nil ) 2 + (z ■ - -U-D 


(7) 


B. Closed-Form Green ’s Functions for Geometries with the Top Ground Plane 

When the top ground plane is present, all of the four coefficient functions Kf are still 
nonoscillatory but contain a pole at y = 0; as a consequence, Kf can no longer be accurately 
approximated with exponential functions [8]. To overcome this difficulty, let us rewrite G in the 
following manner: 

G(y , y|r 0 ) = R m , n G h (/, y|r 0 ) + G (y, y|r 0 ) (8) 

where G h is the spectral-domain Green’s function for a homogeneous medium, i.e., all dielectric 
layers are replaced by the source layer. R m n is a constant which is determined such that G 

contains a pole at y = 0 and G is a well-behaved function without any poles. R m n can be 
obtained either numerically or analytically by taking limits of G and G h as y— >0. Now the 
technique used in the previous section can be applied to obtain the closed-form expression for G 
in the space domain, and the space-domain expressions of G are obtained once the corresponding 
expressions of G ^ are determined. 

Expressions of G h in the space domain can be easily obtained using the image theory approach 
and are given by 
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2q 

Unfortunately, all expressions are written in terms of infinite series. Although G can be 

alternatively expressed using a closed-form formula [8], such closed-form formulas cannot be been 

found for G 3D ' h and G p ' h . Since the closed-form formula for G 2DM requires numerical 

2 D h 

integration when the moment matrix is computed, we will simply use (10a) to evaluate G 

Therefore, an infinite-series expression, in general, cannot be avoided for G 2D , G , and G p 

when the top ground plane is present. However, the expressions for G 2D , G 3D , and G p given in 
this paper are still numerically more efficient than the ones obtained from the conventional image 
method since a nested infinite series of the conventional method is reduced to a simple infinite 
series without nesting as shown in the above equations. 1 For this reason we shall still refer to 
G 2D , G 3D , and G p given by (9), (10a), (10b), and (10c) as closed-form Green’s functions. 


IV. FORMULATION OF AN INTEGRAL EQUATION 


In this section, an integral equation is formulated in terms of the excess charge distribution 
using the closed-form Green’s functions derived in the previous section. Figure 2 shows the planar 
view of the general geometry of a discontinuity, which consists of traces and a junction region. 
This general geometry represents most of the common strip discontinuities, e.g., an open-end, a 
nonothogonal bend, and various junctions. Although the present approach can handle conductors 
with finite thicknesses, the conductor thicknesses are assumed to be infinitely thin in this paper. 
The discontinuity structure is embedded in a layered dielectric medium, which is shown in Fig. 1 . 

1 Alternatively, infinite series can be avoided by modeling the top ground plane as an additional conductor and 

using the Green’s functions in the previous section. 
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The integral equation relating the electrostatic potential 0 and the charge density q on the 
surface of a conductor is given by 

0(r) = J G W (r\r' )q(r' )dr> = (c W ,q) ( 10) 

Q 

where Q are the surfaces of a conductor: traces and a junction region. To simplify the notation the 
integration is symbolically written as ( v >. Now let us rewrite the charge density q in the following 

manner: 

f qj (r), if r is on the junction region 
q(r) = { . OD 

[q l T (r), if r is on the ith trace 

Then, (10) becomes 

0(r) = (G 3D , i; ) = {c 3D ,„) + i(G 3D ,^) (12) 

1 = 1 

Decomposing the charge densities q l j into the uniform charge density qj and the excess charge 
density q e j Cess ' 1 for each trace: 

q‘ T (r) = qf f ’ i (r) + q e T XCeSS ' i (r) (13) 

Here, the uniform charge density q ™ if ' 1 is obtained by solving a 2-D problem, in which it is 
assumed that only the ith trace is present in the medium and that the ith trace is infinitely long in 
both directions. A detailed discussion for solving 2-D problems is given in [8]. The uniform 
charge density exists only on the ith trace, which is a semi-infinite line; hence, G semi 

should be used to compute the potential due to qf#*. Using (13), (12) can be written as follows: 

- (g 30 .9j) + X(g 3D .9t”“") 04) 

i = 1 '=1 

The integral equation (14) can now be solved by using the method of moments. The collocation 
method is used in this paper. The closed-form formula for the integration involving G 30 is given 
discussed in [8], whereas the integration involving G semi can be analytically integrated using the 
following formula: 
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(15) 


Now once ( 1 5) is solved, the excess (equivalent) capacitance c e can be obtained by 


c 


e 



(16) 


where £?/ is the total charge on the junction region, and Qj XCess ' 1 is the total excess charge on the 
ith trace. Throughout the formulation, we have assumed that the junction region exists between 
traces. The formulation for cases without the junction region, such as an open end and step 
junctions, can be easily obtained simply by removing terms corresponding to qj. 


IV. NUMERICAL EXAMPLES 


A computer program is written based on the technique discussed in the previous sections. The 
program assumes that the shape of the junction region, (see Fig. 2), if it exists, is polygonal, and it 
can handle an arbitrary number of dielectric layers as well as traces. 

Excess capacitances for four common strip discontinuities, an open end, a step junction, a 
bend, and a T junction, are computed using the program. The geometries of discontinuities with 
their corresponding equivalent circuits are shown in Fig. 3. The following parameters are used: 1) 
an open end: w = 0.5 mm, 2) a step junction: w, = 0.1 mm and w 2 =0.2 mm, 3) a right-angle 
bend: w, = vv 2 =0.15 mm, and 4) a T junction: w, = w 2 = w 3 =0.15 mm. Three different types 

of media are considered for each discontinuity with the following parameters (see Fig. 5): 1) an 
open end: £, = 4.2, £ 2 = 2.5, y, = 1.0 mm, y 2 = 1.5 mm, and y 3 = 2.0 mm, 2) a step junction: 
e,=6.0, £ 2 = 4.2, y, =0.1 mm, y 2 =0.2 mm, and y 3 = 0.3 mm, 3) a bend: £,=2.5, 
£ 2 =4.2, y, =0.15 mm, y 2 =0.3 mm, and y 3 =0.5 mm, 4)aT junction: £, =2.5, £ 2 =4.2, 
y, = 0.15 mm, y 2 = 0.3 mm, and y 3 = 0.5 mm. All discontinuities are assumed to be embedded 
at y _ y j 0 pi ace 3.0 unknowns for the excess charge distribution, the length of each trace is 
truncated at / = 8w . The total numbers of unknowns per each trace were 50 for a 2-D problem and 
160 for a 3-D problem, whereas 100 unknowns were used for the junction region. The maximum 
number of exponentials used to approximate each coefficient function AT* was 5. 
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The computed results are shown in Table I with the comparison data for a microstrip case (Fig. 
4(a)). A good agreement was found overall as shown in the table. It is interesting to note that for 
some cases the value of an excess capacitance turns out to be negative. Although a physical 
capacitance must be positive, an excess (equivalent) capacitance is hypothetical and can be 
negative. 

V. CONCLUSIONS AND FUTURE WORK 


An efficient method to compute excess capacitances of strip discontinuities was discussed in 
this paper. Complete expressions of closed-form Green’s functions for a point charge, a line 
charge, and a semi-infinite line charge have been derived. Unlike other approaches, only one 
integral equation is employed in this paper to handle various strip discontinuities instead of 
formulating an integral equation for each discontinuity type. The numerical results for a microstrip 
case agreed well with other published results. 
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Fig. 2. General geometry of a strip discontinuity. 
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Table I. The numerical results (Units are in femtofarad). 



Medium 1 

Medium 2 

Medium 3 

Computation 

Others 

Computation 

Computation 

Open End 

17.33 

17.0 [1] 

23.52 

19.62 

Step Junction 

1.120 

1.05 [1], 0.74 [7] 

1.352 

0.609 

Bend 

6.210 

6.75 [1], 5.8 [7] 

7.006 

9.184 

T Junction 

1.385 

1.9 [7] 

-4.917 

-0.818 
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Abs tract This paper presents a method for the circuit simulation of uniform 

multiconductor lossy frequency-dependent lines characterized by sampled 
frequency-domain responses. The implementation includes ac, dc and transient 
analyses. The method combines element characteristics which do not require 
introduction of current variables, open-loop characterization which results in the 
simplest transfer functions, novel frequency-domain matrix rational 
approximation method, novel matrix delay separation technique, and matrix 
indirect numerical integration formulas. The method is reliable, accurate and as 
efficient as the simple replacement of interconnects by lumped resistors. 


I. Introduction 

rT^HE PROBLEM of the transmission line simulation gained special importance with the 
J_ development of high-speed digital electronics. As transient times become faster, the 
transmission line behavior of electronic interconnects starts to significantly affect transient 
waveforms, and accurate modeling of on-board and even on-chip interconnects becomes an 
essential part of the design process. The complexity of contemporary digital circuits calls for 
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accurate transmission line models efficient enough to perform the simultaneous circuit simulation 
of thousands of lossy coupled frequency-dependent lines surrounded by thousands of nonlinear 
active devices. Lines to be simulated may be characterized by measured or electromagneticly 
simulated samples of their responses. 

A substantial amount of study has been devoted to the transient simulation of transmission 
lines in recent years [1]-[19]. This paper presents a method for circuit simulation of uniform 
multiconductor lossy frequency-dependent transmission lines characterized by sampled frequency- 
domain responses. The method is based on the approach described in [20]. The method supports 
variable time stepping and has linear computational complexity. The method has been adopted in 
several industrial and commercial circuit simulators, and, in numerous real-life simulation 
exercises, proved to be reliable, accurate and as efficient as the simple replacement of interconnects 
with lumped resistors. 

The method employs element characteristics which do not require the introduction of 
current variables. As a result, the method increases neither the number of nodes nor the number of 
circuit variables, and the lines are incorporated into the circuit simulation without any increase in 
the circuit solution time. Element characteristics for the ac and dc analyses are based on the Y 
parameters. 

For the transient analysis, the method uses open-loop transfer function characterization 
(direct characterization in terms of the propagation functions and characteristic admittances), which 
opens the feedback loop formed by the reflections from the terminations, and results in the simplest 
characteristic responses. It is shown that the open-loop characterization is equivalent to the 
generalized method characteristic . 

The transient analysis is carried out by matrix indirect numerical integration. It has linear 
computational complexity, and ideal convergence, accuracy and stability properties. 

To apply numerical integration to the propagation functions and characteristic admittances 
given as a set of frequency-domain samples, the method employs difference approximations. It 
fits the samples with a rational polynomial function and expresses the numerical integration 
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formulas in terms of the approximation parameters. Because of the simplicity and aperiodicity of 
the open-loop transfer functions, excellent approximation accuracy is attained with low-order real- 

pole approximation and a few original samples. 

The approximation is performed by novel matrix direct constrained method. The method 
fits arbitrary number of arbitrary-spaced complex samples of a matrix transfer function with the 
matrix rational polynomial functions with poles constrained to the left half of the complex plane to 
insure stability. The method employs no iterative or relaxation techniques, only direct linear 
solutions and a polynomial factoring, and is computationally efficient and free from convergence 

problems. 

Before the approximation, the delay is separated from the matrix propagation functions by a 
novel matrix delay separation technique. The technique eliminates the problems associated with the 
delay separation based on the diagonalization with the frequency-dependent modal transformation 
matrices, which are nonminimum-phase functions of frequency with unstable time-domain 

responses. 

Throughout the paper, capital boldface, small boldface and normal italic symbols will be 
used to denote matrices, vectors and scalars, respectively. Since only multiconductor lines will be 
considered, the modifier “matrix” for the transfer functions will be omitted in the future for brevity. 

II. Frequency-Domain Line Model for Transient 

Analysis 

The frequency-domain element characteristic (for the transient analysis) which does not require the 

introduction of current variables and is suitable for the line modeling is given by 

fi^co) = Y,(co) v,(co) - j,(co) (1) 

[i 2 (a>) = Y 2 (o>) v 2 (cd) - j 2 (to). 

The conventions for the terminal voltages and currents are shown in Fig. 1. The expressions 
relating the matrix admittances Y, and Y 2 and vector current sources j, and j 2 to the 

transmission line characteristics are derived directly from the continuity conditions for the voltages 
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and currents at the line terminals. To separate forward and backward waves and open the feedback 
loop, the current source j, must depend only on the backward wave, and j 2 only on the forward 

wave. This condition uniquely defines Y, , Y 2 and j, , j 2 as follows: 

Y,(co) = Y 2 (cd) = Y c (co) (2) 

and 

jj,((0) = 2i bl ((0) {3) 

jj 2 (co) = 2i f2 (co), 

where Y c stands for the characteristic admittance, and the forward and backward current waves, 
i r , , i f2 and i bl , i b2 are related as: 

Ji bI (co) = W lb ((o) [i b2 (to) = i 2 (o)) + i r2 (co)] (4) 

|i f2 (( 0 ) = W lf (to) [i„ (co) = i,(co) + i b ,(co)]. 

The propagation functions for the forward and backward current waves are equal, 
w ir = w ib = W,. The open-loop device model (l)-(4) is equivalent to the generalized method of 

characteristics [14], [11], [12]. 

The propagation function and characteristic admittance can be computed from the insertion 
loss data [14], scattering parameters [19], or distributed RLGC parameters: 

W,(co) = e' K,< “ M 5 

Y c (co) = K,(co) Z _I (co), 

where l is the length of the line, 

K,(co) = (Y(co)Z(co))i 

is the propagation constant for current waves, 

Y(co) = G(co) + ;coC(co) > 

Z(co) = R(co) + ;'coL(co) 

are the admittance and impedance per unit length, and R, L, G and C are the resistance, 
inductance, conductance and capacitance per unit length. Boldface (.) 2 and e u denote matrix 
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square root and matrix exponential, respectively. 


III. Difference Approximation 

To perform the transient analysis, indirect numerical integration is applied to the propagation 
functions and characteristic admittances in the frequency -domain line model (l)-(4) by using the 
difference approximation method. 

For the difference approximation in the parallel canonic form, samples of the frequency- 
domain transfer function are approximated with the rational polynomial function 


M 

HU©) = HL + X 


f 0 l + ;'CO/(O c 


(5) 


or samples of the time-domain unit-step response are approximated with the exponential series 

h(0 = H 0 -Xc' u '"'A m - 

m = 1 

where H 0 and H„ denote the initial and final values of the approximating transfer function 


H(;a». 

Once the approximation has been performed, indirect numerical integration formulas 
(discrete-time difference equations) are readily given in terms of the approximation parameters. 
For the step invariance the formulas are 

( M 

y(0 = H.x(0 + X z m(0 (6 n 

1.(1, ) = (l - e'“" T - )A. *(»..,) + e '“" T ' z - 

where x, y and z. stand for the excitation, response and state variables, respectively, and 
x = f — f is the time step at the nth transient iteration. 

For the ramp invariance 


y(rJ = H 0 x(r„)-Xz m (r n ) 

m~\ 


( 7 ) 
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where 


D m (rj = 


CD, T 


A„. 


An alternative form of the ramp-invariant indirect numerical integration formula has the coefficients 
of the present-time sample of the excitation lumped together 


M 


y(rJ = fH 0 -XD m (7jW n )-X*m^) 

V ms) J m=1 


( 8 ) 


,(rj = (D " r " - D m {T n ))x(r fl _, ) + )• 


This form is especially suitable for discretization of characteristic admittance, because, for 
admittances, the present- and past-time terms of the numerical integration formulas have different 

physical meaning. 

Before the approximation, the delay is separated from the matrix propagation function 
using the matrix delay separation technique. It represents the matrix propagation function in the 


following form: 

W,(cd) = W,(co)M i e-'“ Tta M;\ 

where T Im and M, are the constant eigenvalue and eigenvector matrices of the propagation delay 
matrix 

T, = (C(°°) L(~))z / = K,(°o) / 


Difference approximation is applied to the delayless propagation function 

W,(to) = W,<ct>) e'" T ‘ 

The delays in the diagonal modal-delay matrix T lm are modeled using a low-order spline of the 
simulated time points. 

The difference approximation is applied to the delayless propagation function and 
characteristic admittance. For the characteristic admittances in (1), the excitations are the terminal 
voltages and the responses are the terminal currents. For the propagation functions in (4), the 
excitations and responses are the current waves. 
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Since the open-loop functions are aperiodic, they have to be approximated with only real 
poles, — to rm . In addition, the poles have to be negative to assure stability. 

To represent the original functions accurately with the minimum number of samples, the 
variation of the original function from sample to sample should be about the same. The following 
empirical formula for the sampling frequencies was found to provide good results 

to* =w^l-cos^j, k = 0 ,\,...,K. 

The end of the approximation interval, to* , should be chosen so that the original function would 
closely approach its final value. This assures that the resulting indirect numerical integration 
formulas will be accurate in the full frequency and time ranges from zero to infinity. 


IV. Matrix Complex Rational Approximation 
Method for Frequency-Domain Difference 

Approximation 

The method fits samples of an N x N complex matrix transfer function H(co) with the rational 
polynomial function (5) at the set of arbitrary spaced frequencies {0, to,, (0 : «*} ■ The method 

proceeds in three steps. 

First, the real part of the sum of the elements of the original matrix function. 


«!«■» = XXlHWlr 


« = 1 j = 1 


is fit with the real part of the complex rational polynomial function, which is a real rational 


polynomial function of squared frequency 


Re(H z (ja») 


c n + c, to 2 + c,(to 2 ) 2 +... + C, (to 2 )" 
l + P,to 2 + (3 2 (to 2 ) 2 + ... + (3 M (or) M 


(9) 


The following linear system of equations results from matching the real part of the original function 
with (9) at the set of frequencies and premultiplying both sides of each equation with the 


denominator 
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"1 0 . 

O 

O 

0 

c o 


Re{H z {0)) ' 

1 co; 

.. CD 2 ” j -C0?te(// s (0),)) •• 

. -CD^fle(// r (CD,)) 

^ M 


Re{H z ( CO,)) 

1 co 2 * . 

.. (0 2M i -CD 2 fo^CD*)) .. 

. -C0 2M /?<?(#£ (CD*)) 

Pi' 


Re{H z ( co*))_ 




JV 




For interpolation, K = 2M and solving (10) produces a rational polynomial function which 
coincides with the real part of the original function at all of the sampling points. For a set of 
samples larger than 2M + 1, the least square solution of the system (10) can be obtained. 
However, it minimizes the approximation error premultiplied with the denominator, which can lead 
to inaccurate approximation. Better results are achieved with the method of averages [21], which 
partitions the larger number of equations into 2M + 1 subsets in order of the increasing of w . The 
equations within each subset are added up, which makes the system consistent. The method is 

effective in averaging out the noise in measured data. 

After the real part has been approximated, the denominator of (9) is factored yielding the 
squared poles, -cd 2 „,. Consequently, no unstable right-half-plane poles can be produced. 

However, there still can be spurious complex conjugate and purely imaginary poles, which are 
removed. The remaining real negative poles are used to formulate the equations for the partial 
expansion coefficients, A m , of (5). As a result, the order M of (5) is less or equal to that of (9). 

Matching the real and imaginary parts of each element of the original matrix transfer 
function H(0)) with the corresponding parts of elements of (5) at the set of frequencies 
{0, to,, 0) 2 ,..., co*} leads to the following linear systems of equations 
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A = 


'1 

1 

1 

1 

1 

l 

1 +03^ /CO", 

1 + co'/cd- m 

1 

1 

1 

l + or*/co;, 

1 + (0^/(0^ 

0 

-0),/(0 cl 

—CO, /CD,* 

l + cof /( 0 *, 

1+C0'/C0^ 

0 

-CDk/CDc 

3 

3 

* 

3 

1 

l + (o;/(o', 

l + fl£/co-,_ 


x = 


[H.r 

M, 

[A,], 

[AmL 


= b = 


Re I 


Re I 
Im 


[H(0)]„ 

([««■>,>]„) 

fHMj 


Im 


([»«■>„>]„) 


(ID 


For interpolation M = 2K . and both real and imaginary parts of the original transfer function are 
matched exactly at all of the K frequency points and dc. For an arbitrary larger number of points, 


the least square solution of ( 1 1 ) is obtained from 

A t Ax = A T b. 

The total computational complexity of the approximation method is that of N +1 real 
linear solutions and one polynomial factoring. The orders of the polynomial and linear systems 
depend only on the order of the approximation and not on the number of the original function 
samples. Since no iterative or relaxation techniques are involved, the method is free from 
convergence problems. 

Fig. 2 shows an example of the fourth-order approximation of an open-loop transmission- 
line transfer function. As can be observed, although only nine samples of the original function 
were used, the approximation exhibits an excellent match in the full frequency range. In general, 
due to their simplicity, the open-loop functions can be accurately fit with the 3rd-9th-order 

approximation. 


V. Companion Model 

By applying the difference approximation to the propagation function and characteristic admittance, 
the frequency-domain element characteristic (1) is transformed into the following discrete-time 
element characteristic, or companion model 
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i l (^) = Y i ajv 1 (/J-j ,(0 (12) 

.i 2 0 n ) = Y 2 (/Jv 2 (rJ-j 2 (rJ. 

The circuit-diagram interpretation of the companion model is shown in Fig. 3. 

The admittances Y, and Y 2 represent present-time coefficients in the indirect numerical 

integration formulas for the admittances Y, and Y 2 ■ The current sources j 2 and j 2 combine the 
currents j Y and j Yj corresponding to the remaining parts of the numerical integration formulas for 
the admittances, and j, and j 2 which are given by the discretized Eqs. (3) and (4) 

j ,(0 = - jY ,(0 + j .(0 (13) 

j 2 (0 = - 

Eqs. (3) and (4) do not contribute to the admittance part of the companion model because the 
propagation functions contain a delay. 

The Modified Nodal Approach (MNA) stamp corresponding to the companion model (12) 


is: 

KCLat 
node : 

v u - vv V r j v 2J - V 2 , | V, 




1.1 

~ 1 ! ! 1 

V,, 


[i>l 

1 .N 

... <>£ 
•l 

i 

tH 

> 



[i 

Y 

- -£I*,L |XX 1 

V'r 

= 

-fiu 

r«) 

2.1 

~ p 1 ' I j-Iiu 

V u 


~ [i:], ' 

i i Y. 1 

1 JL 9 1 N r „ , 

: 



2.N 

-Ilv.L, 

Vu, 

J 

— TT - 

2' 

i i-fFv.i - -iiv.L lilt*.], 

1 1 ,.i .-i 1 •-* J 

v : . 


-Llil 

1 =1 


(14) 


In the circuit simulator during the transient analysis, the lines are represented by the tables of 
numbers (14) which are recursively updated at each time iteration using numerical integration. The 
left-hand side of the stamp (14) has to be updated only when the value of the time step changes. If 



the step-invariant indirect numerical integration formulas (6) are used, the left-hand side of the 
stamp becomes independent of the time step, and only the right-hand side vector has to be updated. 
Since the terminal currents are not introduced as variables, the values of i, and i 2 in (4) are 

computed from (12). 

VI. Line Model for AC and DC Analyses 

For ac and dc analyses, the complexity of the transfer functions is not important, and the element 

characteristic which does not require the introduction of current variables and is suitable for the 

ac/dc modeling of transmission lines is the K-parameter characteristic 

j i,((o) = Y„(co) v,(co) + Y, 2 (co) v 2 (cd) ( , 5) 

ji,(co) = Y 21 (to) v, (co) + Y 22 (co) v 2 (co). 

The T-parameters are related to the open-loop functions as follows 

Y„(co) = Y 22 (co) = Y c (co) + 2[l - W J (co)] ' Wj (co)Y c (co), ( J6) 

Y 12 (co) = Y 2l (co) = -2[l - Wj(co)]' 1 W,((o)Y c (co), 
where I is the identity matrix. The expressions were derived by eliminating j, and j 2 from Eqs. 

(1).(4) > and transforming them to the form of (15). 

The dc model is merely the ac model at zero frequency. For the limiting case of lines with 

zero distributed conductance, G = 0 , the dc values of the K-parameters are 

Y„(0) = Y 22 (0) = -Y i2 (0) = -Y 21 (0) = ~R ' • 

The MNA stamp corresponding to (15) is 
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J 1 l u 


Y - 1 - 9*13 

-i[v„i, - sfiv-irl _^'_ 

'“■ f^.Ci „ "I-flM, ^ 


V i '": ’i Y i '" : 

2.N 121 [-| |Y ^ 1 l “ 

2 ‘ ~XIYji I* - ^[Y l> L|lf(Y 1 ,]J | ^|[V»] 11 - -t [Y “l*]l5 [Y “k. Vr . 



VII. Initial Conditions for Transient Analysis 

The dc model is used to perform the operating-point (op) analysis before the transient simulation. 
The op solution is then used as initial conditions for the transient analysis. The initial conditions 
for the indirect numerical integration are the dc values of the state variables, which are related to the 
dc value of the excitation, * 0 , as follows: z m (t 0 ) = a m x 0 for the step-invariant case (6), z m (t 0 ) = 0 
for the ramp-invariant case (7), and z m (t 0 ) = -d m (T, )x 0 for the ramp-invariant case (8). The dc 
values of i„ and i b2 , which serve as excitations for the propagation functions in (4), have to be 
expressed in terms of the terminal voltages obtained from the op analysis. Resolving Eqs. (l)-(4) 
leads to 

i t ,(0) = [l-Wj(0)]‘ , [V,(0)v,(0)-W 1 (0)Y,(0)v ! (0)] 

‘i bl (0) = [l-W;(0)l‘ , [Y I (0)v ! (0)-W I (°)Y,(0)v,(° ) ]. 

For the limiting case of G = 0 , the expressions become 

i„(0) = -i„<0) = 4 R- 1 [v, (0) - v j(0)] . 



VIII. Simulation Algorithm 

For an MNA-based simulator, the line simulation algorithm is as follows. 
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1 . Before the transient analysis: 

a) Perform op analysis of the circuit to find the initial conditions for the transient analysis. 
Use the ac/dc model (15)-(16). 

b) For each line in the circuit, perform the difference approximation of each element of the 
propagation function and characteristic admittance matrices. 

2 . At each time iteration: 

Recursively update the line stamps using the indirect numerical integration formulas 
obtained at step 1(b) and companion model ( 1 2)-( 14). 

Since the method introduces neither additional nodes nor current variables, the line 
modeling does not increase at all the circuit solution time. The only additional time is required to 
perform a low-order interpolation once in the beginning of the simulation, and for a low-order 
numerical integration. As shown in the next section, this time is negligibly small compared to the 
circuit solution time. 


IX. Numerical Results 

The method has been adopted in several industrial and commercial circuit simulators, and, in 
numerous real-life simulation exercises, proved to be reliable, accurate and efficient. Table II 
presents relative runtime data for circuits of various types and sizes. As can be observed, even for 
the worst case of a small circuit consisting of lines only, the model is virtually as efficient as the 
simple replacement of interconnects with lumped resistors. The resistive model was chosen for the 
comparison because it represents the limiting case in the simplicity and computational efficiency of 
the interconnect modeling. 

Fig. 4 verifies the method’s accuracy with Spice3e2 lossy multiconductor line model [8]. 
A simple network was chosen as an example to reduce the influence of factors other than the line 
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model on the simulation accuracy. A variable, third- to fifth-order frequency-domain difference 
approximation was applied. As can be observed, the compared waveforms are indistinguishable. 
The runtime for the presented method was three orders of magnitude shorter than that for the Spice 
model, which is, in turn, an order of magnitude faster than segmentation models. 

The presented line model has been implemented in a commercial CAD product DFSignoise 
from Cadence. The following examples show results for real circuits analyzed with than, the 
simulator deployed by Cadence's DFSignoise. The examples compare the suggested line model 
with segmentation models. The suggested model included frequency-dependent loss, whereas the 

segmentation models included no frequency dependence. 

The first example is a long cable driven by an ECL buffer at a high, 400 MHz, speed. The 
nonlinear characteristics of the driver and receiver were described in the I/O Buffer Information 
Standard (IBIS) data format [22], Fig. 5 shows the circuit schematic and the IBIS data, and 
compares the simulation results obtained with the lumped RLGC and pseudo-lumped (distnbuted- 
LC/lumped-/?G) segmentation models and the suggested line model. The number of segments for 
the segmentation models was over 400. The automatic algorithms used to determine the number of 
segments took into account the rise and fall times and the expected attenuation. 

The waveform simulated with the suggested line model shows considerable attenuation due 
to the frequency-dependent dielectric loss and conductor loss with the skin effect. As can be 
observed, the segmentation models without the frequency-dependent loss consistently underpredict 
the attenuation and resulting loss of the noise margin, accurate prediction of which is important at 
these speeds. The lumped segmentation model also shows spurious peaks which are artifacts of 

that method. 

Fig. 6 shows a similar comparison between the pseudo-lumped segmentation model and 
suggested line model for a circuit driven by IBIS CMOS buffers at 200 MHz. This network is 
complicated and contains 45 two-conductor, 8 three-conductor and 5 four-conductor lines. The 
longest transmission line in this case was a two-conductor line of 0.5 m. The longest 
multiconductor line was a four-conductor line of 0.26 m. The largest delay was 7.8 ns. In this 
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case, the two models produced comparable results, since the rise time of the driver was 0.9 ns as 
opposed to the 0.3 ns of the ECL driver. The figure also shows the crosstalk waveform V, . The 
crosstalk waveforms for the suggested and segmentation models were very close. 


X. Conclusions 

This paper presented a method for circuit simulation of uniform multiconductor lossy frequency- 
dependent lines characterized by sampled frequency-domain responses. The method uses element 
characteristics which do not require the introduction of current variables, simple open-loop 
transmission line characterization, matrix indirect numerical integration formulas, novel direct 
matrix rational approximation method, and novel matrix delay separation technique. The complete 
step-by-step implementation of the method was presented, including ac. dc and transient analyses, 
extraction of initial conditions from the op analysis, and MNA stamps. 

The method is compatible with recursive time-domain solvers employed by circuit 
simulators and supports variable time-stepping. The method has been adopted in several industrial 
and commercial circuit simulators, and, in numerous real-life simulation exercises, proved to be 
reliable and accurate. It was shown on an extensive set of runtime data that the method is as 
efficient as simple replacement of interconnects with lumped resistors. 

Acknowledgment 

The authors wish to thank Mr. H. Satsangi for installing the presented line model in Spice, and Dr. 
R. Wang for implementing the model in the circuit simulator at Intel. 

References 

[1] T Dhaene, L. Martens and D. De Zutter, “Transient simulation of arbitrary nonuniform 
interconnection structures characterized by scattering parameters,” IEEE Trans. Circuits 
Syst., vol. 39, pp. 928-937, Nov. 1992. 



16 


Syst. , vol. 39, pp. 928-937, Nov. 1992. 

[2] J. E. Bracken, V. Raghavan and R. A. Rohrer, “Interconnect simulation with asymptotic 
waveform evaluation (AWE),” IEEE Trans. Circuits Syst., vol. 39, pp. 869-878, Nov. 
1992. 

[3] T. Dhaene and D. De Zutter, “Selection of lumped element models for coupled lossy 
transmission lines,” IEEE Trans. Computer-Aided Design, vol. 1 1, pp. 805-815, Jul. 1992. 

[4] K. S. Oh and J. E. Schutt-Aine, “Transient analysis of coupled, tapered transmission lines 
with arbitrary nonlinear terminations,” IEEE Trans. Microwave Theory Tech., vol. 41, pp. 
268-273, Feb. 1993. 

[5] J. E. Schutt-Aine and R. Mittra, “Nonlinear transient analysis of coupled transmission 
lines,” IEEE Trans. Circuits Syst., vol. 36, pp. 959-967, Jul. 1989. 

[6] I. Maio, S. Pignari and F. Canavero, “Influence of the line characterization on the transient 
analysis of nonlinearly loaded lossy transmission lines,” IEEE Trans. Circuits Syst., vol. 
41, pp. 197-209, Mar. 1994. 

[7] F. Romeo and M. Santomauro, “Time-domain simulation of n coupled transmission lines,” 
IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 131-137, Feb. 1987. 

[8] J. S. Roychowdhury and D. O. Pederson, “Efficient transient simulation of lossy 
interconnect,” in Proc. 28th ACM/IEEE Design Automation Conf., 1991, pp. 740-745. 

[9] F.-Y. Chang, “Transient analysis of lossy transmission lines with arbitrary initial potential 
and current distributions,” IEEE Trans. Circuits Syst., vol. 39, pp. 180-198, Mar. 1992. 

[10] F. H. Branin, Jr., “Transient analysis of lossless transmission lines,” Proc. IEEE, vol. 55, 
pp. 2012-2013, Nov. 1967. 

[11] F.-Y. Chang, “Transient simulation of nonuniform coupled lossy transmission lines 
characterized with frequency-dependent parameters. Part I: Waveform relaxation analysis,” 
IEEE Trans. Circuits Syst., vol. 39, pp. 585-603, Aug. 1992. 

[12] F.-Y. Chang, “Transient simulation of frequency-dependent nonuniform coupled lossy 
transmission lines,” IEEE Trans. Components Packaging Manufacturing Technology, vol. 
17, pp. 3-14, Feb. 1994. 



17 


[13] J. E. Bracken, V. Raghavan and R. A. Rohrer, “Extension of the asymptotic waveform 
evaluation technique with the method of characteristics. Tech. Dig. IEEE Int. Conf. 
Computer-Aided Design , 1992, pp. 71-75. 

[14] A. J. Gruodis and C. S. Chang, “Coupled lossy transmission line characterization and 
simulation,” IBM J. Res. Develop., vol. 25, pp. 25-41, Jan. 1981. 

[15] A. Semiyen and A. Dabuleanu, “Fast and accurate switching transient calculations on 
transmission lines with ground return using recursive convolutions, IEEE Trans. Power 
Appar. Syst., vol. PAS-94, pp. 561-569, Mar./Apr. 1975. 

[16] T K. Tang, M. S. Nakhla and R. Griffith, “Analysis of lossy multiconductor transmission 
lines using the Asymptotic Waveform Evaluation technique, IEEE Trans. Microwave 
Theory Tech., vol. 39. pp. 2107-2116. Dec. 1991. 

[17] C. R. Paul, “On uniform multimo'de transmission lines,” IEEE Trans. Microwave Theory' 
Tech., pp. 556-558, Aug. 1973. 

[18] L. T. Pillage and R. A. Rohrer, “Asymptotic waveform evaluation for timing analysis,” 
IEEE Trans. Computer-Aided Design, vol. 9, pp. 352-366, Apr. 1990. 

[19] J. E. Schutt-Aine, “Modeling and simulation of high-speed digital circuit interconnections, 
Ph.D. dissertation, University of Illinois at Urbana, 1988. 

[20] D. B. Kuznetsov and J. E. Schutt-Aine, “Optimal transient simulation transmission lines,” 
accepted for publication in IEEE Trans. Circuits Syst.. 

[21] I. N. Bronshteyn and K. A. Semendyaev, Spravochnik po matematike. Moscow: Nauka, 
1967, p. 578. In Russian. 

[22] D. Duehren, W. Hobbs, A. Muranyi and R. Rosenbaum, "I/O-buffer modeling spec 
simplifies simulation for high-speed systems," Electronic Design News (EDN), p. 65, Mar. 

1995. 



18 


Figure Captions 

Fig. 1 • Conventions for the voltages and currents at the line terminals. 

Fig. 2. An example of the 4th-order complex rational approximation. The original function is 
shown by the thin continuous curves and the approximating function by the thick dashed 

curves. 

Fig. 3. Companion model for a transmission line. 

Fig. 4. (a) the network configuration and (b) comparison of the transient waveforms generated 

using the line model installed in an MNA-based circuit simulator (thick broken curves) 
and Spice3e2 (thin continuous curves). R ] =R S = 50 ft, fl 2 =fi 6 =l kft, Rj=R t = 10 Mft; 
self-inductance L,= 418 nH/m. self-capacitance C =94 pF/m, mutual inductance L m =125 
nH/m, mutual capacitance C m =22 pF/m, R= 15 ft/m, G= 0, 1=0.611 m (all signal 
conductors are the same). 

Fig. 5. (a) the circuit schematic, (b) IBIS model for the ECL buffers, and (c) the transient 

waveforms simulated with the segmentation and suggested line models. The line 
parameters are: L= 502 nH/m, C=67.9 pF/m, *={3.25 ft/m at dc, 87.7 ft/m at 10 
GHz}, G={430 fS/m at dc, 4.3 mS/m at 10 GHz}, /=0.88 m. 

Fig. 6. (a) the circuit block schematic, (b) IBIS model for the CMOS buffers, and (c) the 

transient waveforms simulated with the segmentation and suggested line models. 

Table Captions 


Table I. Relative runtimes. 
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